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Measurements of Turbulent Skin Friction on 
Cylinders in Axial Flow at Subsonic and 
Supersonic Velocities 


DEAN R. CHAPMAN* ann ROBERT H. KESTER* 
Ames Aeronautical Laboratory, NACA 


SUMMARY greatly from each other, thereby indicating need for future ex 
we . — periments at hypersonic Mach Numbers. 

The principal results of a study of various theories for calcu 
lating turbulent skin friction in compressible flow are briefly re 
viewed. These theories, 16 in number, predict widely different NOTATION 
effects of Mach Number on skin friction and, hence, also on heat 
transfer. Such discrepancies emphasize the inadequacy of the 
present turbulent boundary-layer theory and the need for ex 


= average skin-friction coefficient, (force)/|[('/2) 
(p.,, U..”) (wetted area) | 
perimental measurements. , = local skin-friction coefficient, dCr/d(x/L) 
Systematic experiments have been conducted to determine dimensionless heat-transfer coefficient (Stanton 


the magnitude of turbulent skin friction along the cylindrical 
portion of insulated cone-cylinder bodies of revolution having 
overall fineness ratios of 10, 15, and 25. Data were obtained by 
directly measuring forces. Boundary-layer surveys were made 
to determine the correction necessary to apply to the force meas 
urements in order to determine the effective starting position of 
the turbulent flow. Mach Numbers between 0.5 and 3.6 and 
Reynolds Numbers between 4 million and 32 million were investi 
gated. At a Mach Number of 2.0, data were obtained for dif 
ferent pressure distributions (by distorting the flexible-plate 
walls of the wind tunnel) in order to evaluate the effect of a mod 


Number ) 
D cylinder diameter 
D,, De = drag of cone-cylinder and cone, respectively 
Da, Dp base drag of cone-cylinder and cone, respectively 
f friction force on cylinder of length / 
Af = friction force on cylinder of length A/ 
F = friction force on cylinder of length L,(F = f + 

Af) 

length of cylinder over which friction force was 
measured 
length of cylinder necessary to produce bound 
erate pressure gradient on turbulent skin friction. ary layer of momentum thickness 6 at Reyn- 
The results obtained with different pressure distributions olds Number UoA//vo 

showed no appreciable effect on average skin friction of moderate 
changes in pressure gradient. At both subsonic and supersonic 
velocities, the skin-friction coefficient was observed to depend 
only to a small extent on the body fineness ratio. For each fine 
ness ratio, however, the effect of Mach Number was found to 
be large, amounting to approximately a 50 per cent reduction in 
skin-friction coefficient as the Mach Number is increased from 
0to4. This effect of Mach Number did not depend on Reynolds 
Number or body fineness ratio and was found to be in good agree- 


corrected length of cylinder, / + Al 
Mach Number 

static pressure 

Prandtl Number 

temperature recovery factor 
Reynolds Number 

temperature 

velocity 

boundary-layer thickness 


ment with the experimental measurements made on flat plates boundary-layer momentum thickness 


by previous investigators. Several of the various published 
theories of turbulent skin friction agreed with each other and 
with the measurements to within about 5 per cent at all Mach 
Numbers investigated. At a Mach Number of 10, however, 
these same theories predict values of skin friction which differ 


—— 


exponent in viscosity-temperature relationship 
(u ~ T°) 
coefficient of viscosity 
y - kinematic viscosity, u/p 
p mass density 
Presented at the Aerodynamics Session, Twenty-First Annual Subscripts 
Meeting, IAS, New York, January 26-29, 1953 
* Aeronautical Research Scientist. 


i = incompressible flow 
0 = conditions at station 0 (see Fig. 3) 
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Ma 
Theories of turbulent skin friction for compressible flow 
over an insulated flat plate. 


Fic. 1. 


© = conditions at outer edge of boundary layer averaged over 
length L 
w = conditions at wall 


INTRODUCTION 


QO’ THE VARIOUS PHASES of experimental fluid me- 
chanics that are of current interest to aerodynam- 
icists, it appears that experiments on turbulent skin 
friction probably have the most extensive history. 
Such experiments have long been important to naval 
architecture, with the result that measurements of sur- 
face friction date at least as far back as the experiments 
of Beaufoy in 1793. A good review of experiments 
prior to 1932 may be found in a paper by Schoenherr.'! 
Theoretical investigations of turbulent skin friction, on 
the other hand, are of comparatively recent origin. 
It is well known that the semiempirical mixing-length 
theories of von Karman and Prandtl, advanced ap- 
proximately 20 years ago, yield accurate predictions of 
turbulent friction for a flat plate in low-speed flow. 
Interestingly enough, though, subsequent detailed hot- 
wire measurements, such as those of Liepmann and 
Laufer? and Townsend? have indicated the fundamental 
assumptions of mixing-length theory to be incorrect. 
The good agreement between predicted skin friction 
and experiment can be attributed in part to the fact 
that mixing-length theory yields arbitrary constants of 
integration which are free to be adjusted by comparison 
with experiment. 

Unfortunately, when mixing-length theory is applied 
to compressible flow, no additional arbitrary constant 
of integration appears, so that adjustment of constants 
of integration to fit compressible-flow experiments is 
not possible. Also, several uncertainties arise in the 
application of mixing-length theory to compressible 
flow because consideration must be given to new vari- 
ables such as the density mixing length and to the vari- 
ation of fluid properties across the boundary layer. 
Since no unique way of precisely accounting for such 
variables is known, mixing-length theory for com- 
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TABLE | 
Values of Re, w, and r for Curves in Fig. 1 





Refer- Re 
Theory ence r w (Millions 
Clemmow II 7 1 0.76 10 
Van Driest II > ok Ws 
Li-Nagamatsu (a = 1) 16 
Modified Frankl-Voishel 13 I l 7 | 
Rubesin 13 I 1 7 
Clemmow III \ 7 l 0.768 11 
Li-Nagamatsu (@ = !/2) p* 16 
Ferrari f 18 
Smith-Harrop 12 l 0.75 10 
' 
Eckert 8 | 0.8 Arbitrary 
Li-Nagamatsu (a = A*) 16 l 0.768 11 
Van Driest I | 10 
Clemmow I > 7 ] 0.76 10 
Li-Nagamatsu (@ = o)§ 16 
Wilson 9 0.89 0.76 10 
Monaghan 11 1 0.76 Arbitrary 
Cope 6 1 0.76 10 
Extended Frankl-Voishel 5, 13 1 1 7 
Tucker II 15 | I Arbitrary 
Young-Janssen 19 7 
von Karman 1 l 0.76 10 
Tuckér | 14 I l Arbitrary 


* Same theory advanced by different investigators. 


pressible flow has yielded a surprising number of greatly 
different predictions of turbulent skin friction. This 
fact is evident from a glance at Fig. 1, which shows all 
theoretical predictions of skin friction for an insulated 


plate known to the present authors.‘~'® The particu- 


lar values of Reynolds Number Re, recovery factor r, } 


and viscosity-temperature exponent w employed for 
each theory are listed in Table 1. At a Mach Number 
of 5, the skin friction predicted by the various theories 
in Fig. | differs by a factor greater than 3, and the differ- 
ence increases at higher Mach Numbers. Only a small 
portion of this difference can be attributed to the 
somewhat different values of Re, 7, and w employed. 
The inadequacy of existing theories in predicting com- 
pressible turbulent skin friction is evident. 

Apart from the obvious direct application of skin- 
friction data to drag calculations, there is an additional 
application that warrants mention since it is becoming 
namely, the application 
Existing theories of heat 


of prime practical importance 
to heat-transfer calculations. 
transfer in turbulent flow are based on an assumed 
analogy between heat and momentum transport. Ac- 
cordingly, these theories generally do not calculate 
the heat transfer directly but calculate instead the ratio 
of heat transfer to skin friction. This is illustrated in 
Fig. 2, where the ratio 2c,/c,, as predicted by various 
analogy theories,”’~*4 is plotted as a function of the fluid 
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MEASUREMENTS OF TUR 


It is seen that the various theories 


Prandt! Number. 
yield considerably different estimates of the heat-transfer 


rate but that, regardless of which theory is used, a cal- 
culation of heat transfer first requires a knowledge of 
skin friction. * 

In cognizance of the limitations of turbulence theory 
and the necessity of knowing the skin friction before 
either drag or heat transfer can be estimated, numerous 
measurements of skin friction have been made during 
the last few years. In this connection, reference is 
made to the measurements of Wilson;’ Rubesin, May- 
dew, and Varga;'* Liepmann and Dhawan;”> Coles;** * 
Brinich and Diaconis;** Bradfield, DeCoursin, and Blu- 
mer;? Cope;*” and Weiler and Hartwig.*! The flat- 
plate measurements of Coles are especially noteworthy. 


The purpose of the present investigation was to meas- 
ure the average turbulent friction on the external sur- 
face of relatively long cylindrical bodies of revolution 
in both supersonic and subsonic flow. It was con- 
sidered desirable to determine the friction by measur- 
ing forces directly, rather than to deduce the friction 
pitot-pressure 


force indirectly from boundary-layer 


surveys. For long cylinders the boundary-layer thick- 
ness becomes comparable to the cylinder radius, under 
which conditions appreciable departures from flat-plate 
friction values would be expected according to the anal- 
yses of Jakob and Dow” Eckert.** Conse- 
quently, in the present tests various cylinders were in- 


vestigated ranging in length from 8 to 23 diameters. 


and 


The present investigation was initiated in 1949, 
though the experimental method finally employed was 
not developed satisfactorily until a year later. The 


data presented herein were obtained during 1950-1952 
in the NACA Ames 1I- by 3-ft. supersonic wind tunnels 
Nos. | and 2 and are intended to summarize the main 
portion of data. It is planned to publish a detailed re- 


port of the investigation at a later date. 


EXPERIMENTAL METHOD 


The general method employed to measure surface 
friction utilizes two foredrag measurements, as illus- 
trated in Fig. 3. The first drag measurement D, and 
corresponding simultaneous base-drag measurement 
Dz, determine the foredrag (D; — Dz) of a cone- 
cylinder as indicated by the top sketch of Fig. 3. A 
second similar drag measurement, as indicated by the 
sketch second from top in Fig. 3, determines the fore- 
drag of the conical nose (D» — Dye). The difference in 
foredrag represents the friction force f acting over the 


cylindrical portion of length /. Thus, 


f = (D, — Dm) — (D2 — Dr) 


*It is necessary, of course, also to know the Prandtl Number 
Thus, heat-transfer calculations at extremely high Mach Num 
bers are limited also by the uncertainty of the value of the Prandtl 


Number for high temperatures 
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Fic. 3. Sketch illustrating method of determining friction force 
and effective starting position of turbulence. 


As indicated in the center sketch of Fig. 3, the 
boundary-layer thickness at the leading edge of the 
cylinder is not zero but is some value 6). Boundary- 
layer surveys in the plane of station 0, as illustrated by 
the sketch second from bottom in Fig. 3, showed that 
5) depended on J/.., Re, nose shape, ete. In order to 
convert force measurements to conditions correspond- 
ing to zero boundary-layer thickness at the leading edge 
of a cylinder, one can imagine the cylindrical surface 
extended upstream a short distance A/ (see bottom 
sketch in Fig. 3) such that the boundary-layer thickness 
developed over the length A/ just matches the bound- 


ary-layer thickness 6). Since boundary-layer surveys 
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Fic. 4. Photographs of typical models. 
WIRE TRIP 
as 
Fic. 5. Photograph showing boundary-layer trip and surface 


Pitot-pressure probe. 


determined the momentum thickness 4) corresponding 
to 6, the small friction force Af that would act over 
the length A/ is readily determined from the equation 


Af = tDpyl 0° (1) 


which represents conservation of momentum. The 
length A/ is determined from an equivalent form of the 


momentum equation 


Al = 20)/Cr, (2) 


Actually, the skin-friction coefficient Cr,, based on the 
Reynolds Number lA//, is not known until A/ is 
known. Hence, the method of determining Cr, and 
A/ was in principle an iterative one, starting with the 
equation 


Cr, = Cr; (C:/Cr;) (3) 


where Cr;, is evaluated at a Reynolds Number lA/ 


vw from the Karman-Schoenherr equation 


0.242/V Cp, = logw (ReCp;) (4) 


and the ratio (Cr/Cp;) is evaluated at an arbitrary 
Reynolds Number from a plot of uncorrected force 
data (this latter ratio does not depend significantly on 
Reynolds Number). It turned out that the correc- 
tions A/ and Af were sufficiently small, and the correc- 
tion method so rapidly convergent, that iteration was 
For each \/.. and Re investigated, Af 


unnecessary. 
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and A/ were calculated from individual boundary-layer 
surveys determining 4). 





The average skin-friction coefficient Cr and corte. 
sponding Reynolds Number Re employed throughout 
are based on the total length L = / + Al and the totaj | 
force f + Af according to the defining equations 


7 f+ af | 
* (1/2)p. Un? [rD(l + Al)] 
Re = U,(1 + Al)/». 


Reference quantities such as M/.., p., U., and v., cor 


respond to conditions at the outer edge of the boundary 





layer averaged over the length L. 
Since it is not possible in any experiment to obtain | 
fully developed turbulent flow starting precisely at the 
leading edge, some theoretical method of determining 
the effective starting position of turbulence inevitably 
is required. The method described above is based on 
the assumption that the ratio Cr/Cpr,; does not depend 
significantly on Reynolds Number. Other methods 
have been employed which involve different assump- 
tions.® !* * 29 Since none of the methods is exact, 
it is necessary in order to obtain accurate data that the 
starting-length correction be sufficiently small so that 
various correction methods would yield substantially 
the same end result. Such conditions were met in the 
present experiments inasmuch as the starting-length 
coefficient was _ relatively 


correction to skin-friction 


small. This correction generally decreased with in- 
creasing Mach Number, varying between 10 and 3 per 
cent for L/D = 8, 6 and 2 per cent for L/D = 13, and 
tand 1 per cent for L/D = 23. Thus, in the most ex 
treme case, a 10 per cent error in the correction itself 
would introduce | per cent error in Cp. 


Photographs of two models employed are shown in 
Fig. 4, and a photograph showing a boundary-layer 
probe, together with a ring trip, is shown in Fig. 5 
The cylinder diameter for each model was | in., and 
Three different 
lengths were investigated (10, 15, and 25 cal., as meas- 
ured from tip of cone to base of cylinder). In the final 
Cr data, these lengths corresponded to L/D values oi 


the included cone angle was 20°. 





approximately 8, 13, and 23, respectively. 

Preliminary to the main set of measurements, 4 | 
series of surface-probe pressure surveys were made in | 
order to find a satisfactory method of tripping the | 


boundary layer. The results for \/,, = 2.9 and for the | 
lowest tunnel pressure employed at this gree), 
Mach Number are shown in Fig. 6. The Mach Num- 

ber (\/),¢¢, used in the ordinate of this figure is the local 
value that would exist on the surface in an inviscid 
flow, as determined by the method of characteristics 
On the smooth body, transition occurred betweet 
about 4 and 7 in. from the cone tip. Consequently, 
some trip was needed to obtain a completely turbulent 


boundary layer over the cylinder. The three trips 
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shown in Fig. 6 were rings, the cross section of which was 
approximately a semicircle. Each trip protruded about 
0.010 in. (radius of the semicircle) from the surface. 
It is seen that, with this thickness of trip, the * ,-in. 
diameter ring moved transition upstream some but not 


enough. The ' 2-in. diameter ring appears sufficient, 
and the * js-in. diameter ring appears more than suffi- 
cient. 


Although the !/2-in. ring was adequate at J/.. = 
99, it was not adequate at \/, = The */;6-in. 
ring proved to be an adequate trip for every Mach 
Number and Reynolds Number investigated in the 


3.6. 


No. 2 wind tunnel; hence, it was selected for use in all 
measure:nents of turbulent skin friction conducted in 
that particular tunnel. In the No. 1 wind tunnel, 
where the turbulence level is much lower and where 
lower Reynolds Numbers are obtained, it was found 


! 


necessary to employ two rings of diameter, * j, and | ; 


in., each protruding 0.005 in. from the surface. 
RESULTS AND DISCUSSION 


Measurements Testing Validity of Experimental Method 


During the course of the main investigation, it was 
believed desirable to conduct a number of supplemen- 
tary experiments. These supplementary tests are dis- 
cussed prior to the presentation of the main body of 
results because they provide some insight into the gen- 
eral validity of the final data obtained. Some of the 
supplementary tests were the following: At \/,. = 
0.81, an ogival nose with rounded tip was investigated, 
in addition to the conical nose, to see whether the data 
depended significantly on nose shape; at J/.. = 2.0, 
one model was tested without a trip, as well as with 
various types of trip, to see whether the data depended 
appreciably on the use of a trip or on the type used; 
also, at /,, = 2.0, tests were made in two different 
wind tunnels (the No. 2 tunnel being of the blowdown 
type with relatively high turbulence level and the No. | 
tunnel being of the continuous-operation type with 
relatively low turbulence level) in order to see whether 
the data were affected by the particular wind tunnel. 
It will suffice here to say that no significant effect on 
Cy of nose shape, type of trip employed, or wind tunnel 
was tound. 

One supplementary investigation, not listed above, 
was made that warrants some discussion. At M/.. = 
2.0, the flexible-plate nozzle walls were purposely dis- 
torted in order to obtain friction data with various 
pressure distributions along the cylinder. 
distributions, plotted in the form of Mach Number dis- 
tribution, are shown in Fig. 7. The corresponding skin- 
friction data are also shown. Evidently, the average 
turbulent skin-friction coefficient is not significantly 
affected by variations in pressure of the amounts inves- 
tigated. It is emphasized that the two pressure dis- 
tributions shown differ from each other by an amount 
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skin-friction coefficient 


much greater than the difference between the normal 
pressure distributions of the present investigation and 
a constant pressure. Also, 1t may be noted that the 
two pressure distributions shown differ from each other 
roughly by the same amount that the pressure dis 
tribution on a 3 per cent thick biconvex airfoil differs 
from a constant pressure; or, alternately, roughly by 
the same amount that the pressure distribution on a 
parabolic-are body of revolution of fineness ratio 15 
differs from a constant pressure. The range of pres- 
sure variation covered by the present tests, therefore, is 
of some practical significance. Since no effect on Cp of 
variations in pressure distribution was observed, the 
skin-friction data to be presented can be considered 
directly applicable not only to the ideal case of con- 
stant pressure but also to slender airfoils. (These data 
can also be applied to pointed bodies of revolution; 
however, strictly speaking, small corrections should be 
made to allow for the relative thinning of the boundary 
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layer over the expanding nose and for the departure 
from two-dimensionality over the aftersection of the 
body due to the boundary-layer thickness becoming 
comparable to the body radius.) 


Measurements of Skin Friction for Various Values of L/D, 

M., and Re 

A plot showing skin-friction coefficient as a function 
of Reynolds Number for the three cylinders (L/D = 
8, 13, and 23) at Mach Numbers of 0.81, 2.5, and 3.6 is 
presented in Fig. 8. At several values of Re the curves 
for different fineness ratios overlap, indicating an in- 
crease in skin’ friction with increasing length-diameter 
ratio. This effect is small, however, amounting to 
slightly less than 5 per cent for L/D = 
with L/D = 8. 


23 compared 


It is noted that in Fig. 8, as well as in all other figures, 
the skin-friction coefficient Cr; for zero Mach Number 
arbitrarily is defined as the value given by the Karman- 
Schoenherr equation [Eq. (4)].. The various curves 
in Fig. 8 appear to be nearly pacallel to the curve for 
M.. = 0. Since a logarithm scale is employed, it fol- 
lows that the ratio Cr/Cp, is approximately independent 


Cr 


Fic. 10. 
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of Reynolds Number. This result can be seen more 
clearly from Fig. 9, where the ratio Cr/Cp; is plotted 
versus Reynolds Number. In this figure, tagged sym. 
bols denote measurements in the No. 1 wind tunnel, 
and untagged symbols denote measurements in the No, 
2 tunnel. No substantial effect of Re on Cr/ Cp, is ob- 
served within the Reynolds Number range investigated, 
Within this range, therefore, it is a good approximation 
to regard Cr/Cp; for a given body as a function of V/ 
alone, provided, of course, that both Cr and Cp, are 
If the Reynolds Number 
range were increased many times, a small effect of Re 


evaluated at the same Re. 


on Cpr/Cpr; might be expected. 

A plot of Cr/Cr; against /., is shown in Fig. 10 
Again, the data generally exhibit a small effect of in- 
creasing skin friction with increasing length-diameter 
tatio. With regard to theoretical calculations of the 
effect of L/D on skin friction, it is to be noted that the 
analyses of Jakob and Dow* and Eckert**® predict an 
increase of skin friction with increasing L/D ratio. 
The analysis of Jakob and Dow, which assumes that 
the rate of boundary-layer growth along a cylinder is 
identical to that along a plate, predicts approximately 
15 per cent increase in friction for L/D = 23 as com- 
pared with L/D = 8. 
assumes that the relation between friction coefficient 


The analysis of Eckert, which 


and boundary-layer thickness is identical for cylinder 
and plate, predicts only 1.5 per cent increase. Al- 
though neither analysis appears quantitatively accu- 
rate, the computed values bracket the corresponding 
experimental value of approximately 5 per cent. 

Thus far, the reference value employed for Cp; is the 
Karman-Schoenherr value for flat-plate flow. Actually, 
the correct reference value for Cr; would depend some- 
what on the L/D ratio of the cylinder. Consequently, 
the true effect of 1/7., on Cp would have to be evaluated 
by using values of Cp, appropriate for each value of 
L/D. of skin-friction 
coetlicient for each cylinder in incompressible flow was 


An experimental value (Cpr,)exp 


determined from the two subsonic measurements by 
applying a 2.5 and 6.1 per cent correction to the data 


1.0 








Mo 


Effect of Mach Number on Cr/Cp; for various length- 
diameter ratios. 
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MEASUREMENTS OF TURBULENT SKIN 


for M.. = 0.51 and 7. = 0.81, respectively, and then 
averaging the two values for each cylinder. These 
corrections, which were applied in order to allow for the 
small effect of compressibility at subsonic velocities, 
were determined by averaging the predicted values 
of Cr/Cr; given by the equations of Cope, Tucker IT, 
and Wilson. Values of (Cp,)¢x». deduced in this manner 
for L/D ratios of 8, 13, and 23 are 2.5, 0.2, and —0.6 
per cent lower, respectively, than the value predicted 
by the Karman-Schoenherr equation. With (Cpr,)exp 
so determined, a plot of Cr/(Crjexp. versus M.. was 
made. This plot showed little scatter (mean deviation 
0.7 per cent) and did not show any discernible differ- 
ence between the various L/D ratios. Consequently, 
the data for the three L/D ratios were averaged at each 
Mach Number. The resulting values, tabulated in 
Table 2, can be compared directly with flat-plate data. 

The present data are compared with data froni other 
experiments in Fig. 11.* The Reynolds Numbers of 
the various experiments are listed in Table 3. Except 


* The measurements of Bradfield, DeCoursin, and Blumer?’ on 
cones are not included in Fig. 11. It is noted, however, that 
these data for cones indicate higher friction than the data in Fig. 
11 for plates and cylinders. 

Measurements of friction on wind-tunnel walls (Cope,* 
Weiler, and Hartwig*'), for which the length Reynolds Number 
was not determined experimentally, are not included in Fig. 11. 

In the experiments of Coles,?* the effective starting position of 
turbulence was not determined. The starting-length correction 
probably would not exceed a few per cent, since the boundary 
layer was tripped near the plate leading edge and the data shown 
are for the station farthest downstream (24 in. from leading edge). 

In the report of Brinich and Diaconis,” the final data were 
plotted versus Reg rather than Re,. In order to convert these data 
to Rey, the effective start of turbulence for each cylinder at each 
tunnel pressure was determined from their measurements of @ at 
upstream stations by the method described in this paper. Their 
measurements of @ at downstream stations, for which the starting- 
length correction to Reynolds Number was less than about 20 
per cent, were then employed to compute Cr. A total of 14 val- 
ues were determined in this manner, the average of which is 


plotted in Fig. 11 


TABLE 2 
Values of Cr/( Cri exp 


M Cr/(CFi)exp 
0.51 0.985 
O.81 0.929 
1.99 0.746 
2.49 0.671 
2.95 0.623 
3.36 0.578 
3.60 0.551 


TABLE 3 
Values of Re for Data of Fig. 11 


Re X 1076 
1 (Approximate) 


Experiment 
Liepmann-Dhawan 
Coles 


Chapman-Kester 6 to 16 
Wilson 10 
Rubesin-Maydew-Varga “§ 
Brinich- Diaconis 3 to 18 
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for the present measurements and the Liepmann- 
Dhawan measurements wherein experimental reference 
values of skin friction for M.. = 0 were deduced, all 
data points shown in Fig. 11 have been referred to the 
incompressible values given by the Karman-Schoen- 
herr equation.t Filled symbols denote data obtained 
by surveying the boundary layer and then calculating 
the friction coefficient by the usual momentum method. 
Unfilled symbols denote data obtained by directly 
The two methods of determining 
It appears 


measuring forces. 
friction yield somewhat different results. 
that the momentum method at certain Mach Numbers 
can yield values on the order of 5 per cent greater than 
corresponding values obtained by the method of di- 
rectly measuring forces. There are at least two possible 
reasons for this: First, the reading of a Pitot tube near 
the surface where the maximum momentum decrement 
occurs is not entirely certain because of the probe inter- 
ference; second, the conventional momentum method 
neglects the contribution of the Reynolds normal stress 
to the momentum integral. 

Three theories that agree with the various measure 
ments within the limit of possible experimental uncer- 
tainty are included in Fig. 11. In making numerical 
computations, the corresponding three authors selected 
different values of w, 7, and Re, as indicated in Table 1. 
For consistency, these three theories have been re- 


computed using the values Re = 10’, w = 0.76, and 
r = 1. Similarly, all three were recomputed for Cpr 
given by the Karman-Schoenherr equation. As a re- 


sult of the recomputation, the curves in Fig. 11 differ 
somewhat from the corresponding curves in Fig. 1. 
t may be noted that the ‘‘extended Frankl-Voishel” 
theory'® also would agree well with experiment up to 
+ The value of c;; corresponding to the Karman-Schoenherr 
equation is 
cy; = dCr;/d(x/l) = 0.557 Cr; /(0.557 + 2V/ Cr;) 


where Cr; is defined by Eq. (4) 
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the maximum Mach Number for which data are now 
available. The extended Frankl-Voishel theory, how- 
ever, is not shown in Fig. 11 because of the time-con- 
suming calculations that would be required to recom- 
pute Cy for the above-mentioned conditions. Cope and 
Tucker each have advanced two theories. The par- 
ticular curves shown in Fig. 11 represent Cope’s ‘‘log- 
law’’ theory and Tucker's theory wherein the mean 
temperature (7 + 7..)/2 is employed as reference 
temperature. 

Although the three theoretical curves shown agree 
fairly well up to \/.. = 4.5, they begin to diverge at 
higher Mach Numbers. In view of this, the upper 
limit for which Cy can at present be estimated to within, 
Clearly, additional 
experiments are needed before turbulent skin friction at 


say, 5 per cent, is about /@, = 5. 


Mach Numbers above this value can be predicted with 


confidence. 
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Structural Analysis and Influence Coefficients 


for Delta Wings 


SAMUEL LEVY? 


National Bureau of Standards 


SUMMARY 


A method is presented for determining the stresses and de 
formations in delta wings of multispar construction having chord 
wise ribs and thin cover sheets. The spars need not be parallel, 
and the cover sheets can have large cutouts. The individual 
spars may be clamped, simply supported, or free at the root 
The method is based on the interaction between the bending 
stiffaesses of the spars and ribs with effective cover sheet and 
is well 


the torsional stiffness of the cover sheet. The method 


suited to the use of high-speed computing machines 


INTRODUCTION 


Ww THE DEMAND FOR increased speed in aircraft, 
aerodynamic considerations required first a 
change in wing shape from the familiar thick, high as 
pect ratio, straight wing to the thin, high aspect ratio, 
sweptback wing. Continued need for even higher 
speeds has led to further modification of the wing con- 
tour, and, as a result, the delta wing with its low aspect 
ratio is coming into wider use. For swept and delta 
wings, conventional methods of analysis using beam 
theory do not give satisfactory results. In the past 
few years, newer types of analysis based on energy 
methods, on plate theory, on continuity considerations, 
and on the assumption of constant chordwise slope of 
multispar wings have been developed.'"> These 
methods have been found generally to give satisfactory 
results for sweptback wings. In some cases, however, 
particularly for multispar wings having delta or other 
low aspect ratio plan forms, a more refined method is 
desirable for a final stress analysis and for the detail 


Asa 


result, consideration has been given at several laborato- 


analysis of ribs and of shear-transmitting fittings. 


ries to methods of analysis which combine properly the 
interaction of ribs, spars, and torque cells. The method 
presented herein is of the latter type. It is intended to 
give improved accuracy and, in addition, to take ac 
count of secondary effects such as chordwise bending, 
discontinuities in the cover sheets, and the actual sup- 
port condition at the root of the individual spars. In 
the case of wings having a wing-fold joint, the last ef- 
fect is particularly important, since the main spar m7y 
be clamped at the fold, while other spars may be simply 
supported or, in some cases, free. 

Received November 10, 1952. 
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The method will be presented in a form particularly 
suitable for use with automatic digital computing equip- 
ment, since such equipment is finding increased use by 
felt that the 
equipment is imperative if the analysis time for highly 


aircraft companies. It is use of such 
redundant structures, such as delta wings, is to be kept 
low enough to allow several alternative configurations 
to be analyzed in a reasonable length of time. 


STRESS IN COVER SHEET 


It is common practice in the analysis of monocoque 
structures having thin cover sheets to consider the 
axial-load and shear carrying capacities of the cover 
sheet separately. Thus, the axial load carrying stiff 
ness of the cover sheet of a conventional wing is ade- 
quately considered by computing // or, in more re- 
fined analyses, by lumping the various portions of the 
cross-sectional area of the cover sheets with that of the 
neighboring spar caps and other axial reinforcements. 
The shear carrying stiffness of the cover sheet is con- 
sidered by computing G/ or equivalent torsional stiff- 
ness expressions. This practice is based on the implicit 
assumption that the stress in the cover sheet in the di- 
rection perpendicular to the spar caps and axial rein- 
forcements is negligible. For high aspect ratio wings 
having relatively little taper or sweepback and having 
bulkheads or ribs perpendicular to the beam axis, this 
implicit assumption is nearly satisfied, and thus theories 
based on it give results that agree well with experiment 

In the root portion of sweptback wings, the assump- 
tion that the stress transverse to the spars and axial 
reinforcements can be neglected must be re-examined. 
For sweptback wings, the stress in the root region can 
be considered to be negligible either (1) normal to the 
spars or (2) parallel to the fuselage. 
(1) in an analysis since it applies also to the outboard 


It is easier to use 


portions of the wing and has been found to give good re- 
sults. ! 

In a multispar delta wing or other wing of low aspect 
ratio, the selection of a direction of negligible stress in 
the cover sheet is no longer possible, and this attractive 
assumption cannot be made without the likelihood of 
substantial error. The angle that the individual spars 
make with the fuselage will probably vary from the 
front to the trailing spar, and the fuselage and ribs may 
not be at right angles to the spars. It is therefore de- 
sirable to consider the load-carrying ability of the cover 
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increase of about 10 per cent for ordinary values 9 
Poisson's ratio, v. If the ribs are extremely light, y taken 
stresses will result, but the curvature in the transverg, proper 


direction will be (— v) times that in the spanwise dire} !ecUV® 
tion, about —0.3 for ordinary values of Poisson's side of 
ratio, v. For intermediate values of rib stiffness, both the - 
of these effects will be smaller. It is felt that, in com. releren 
parison to other uncertainties in the analysis of delt, width 1 
wings, the effect of Poisson's ratio can be omitted for the} oon 
present. choice | 
> In 

METHOD loads 

spars, ) 


The method of analysis is one of consistent deforma points 
tions. The load carried by the wing is considered to be | tance 
the sum of (1) that carried by the spars and cover sheet | with t 
in spanwise bending, (2) that carried by the ribs and] width 
cover sheet in chordwise bending, and (3) that carried | sented 
by the cover sheet in torsion. The load carried by each | sheet ¢ 
of these three mechanisms is determined in terms of the has th 





: ; deflections at the junctions of ribs and spars. The} mum : 
Fic. 1. Lumping of cover sheet area with spars. Dotted i ; : : — , ; 
lines fall midway between spars. The bending stiffness of each Sum Of these loads, expressed in terms of deflections, is} gdere« 
spar is computed assuming the sheet between dotted lines acts equated to the applied external forces at the junctions This n 
with the included spar. Jan ‘ = ‘ j 
rhe equations that result are solved for the deflections }  jnyolv 
at the junctions. From the deflections, the bending} felt th 


sheet both in the direction parallel to the spars and in a ; : ; 
and torsional stresses are in turn determined, thus giv-]  gerody 


transverse direction. 
Because of the effect of Poisson’s ratio, stresses will 
arise in the cover sheets from the restraint of transverse 


ing a complete deformation and stress analysis of the} matior 
wing. asa W 
curvature by the ribs. If the ribs are extremely heavy The ability of the thin cover sheets to carry spanwise | lack 0 
and closely spaced, they will entirely prevent lateral bending load is considered by lumping the sheet with] the ot! 
expansion of the sheet, with the result that the apparent 
stiffness of the sheet in the direction of the spars will be 


the individual spar caps as shown in Fig. 1. The bend-] small 
ing stiffness of each spar is computed assuming the] inters 
sheet to the mid-spar distance to either side as fully ef. 





increased to 1/(1 — v*) of the value with no restraint,an °"** Wh 
fective. 
tion n 
The ability of the thin cover sheets to carry chord} ji: 
wise bending load is considered by lumping the sheet} jiych 
with the individual rib flanges as shown in Fig. 2. The} asin 
bending stiffness is computed assuming full effectiveness wing 
for the sheet up to 0.181 of the rib length to either side deont 
of the rib (or midway between ribs in those cases where | 4, ¢h, 
ribs are spaced more closely than 0.362 of the mb] 
length). The coefficient 0.181 was selected as a reason | 5.44; 
able value for ribs having no bending moments at their 
ends and a gradual increase in bending moments to a Tor 


maximum value near their mid-lengths, on the basis of | torqu 
effective width formulas for wide flange beams (refer- | Jecen! 
ence 9, p. 161). This method of considering the ability | ‘orqu 
of the sheet to carry chordwise bending load is ap- | Stifin 
proximate both in the factor, 0.181, and in the fact that | Core 
the ribs may not be at right angles to the spars, thus | ‘rsio 
violating the conditions of orthogonality between the } Sume 
axial stresses in the sheet resulting from rib and spat This ; 
bending. These errors, though significant insofar as torqu 
the ribs are concerned, are considered to be minor for | ™ the 





the wing as a whole because the loads are not ordinarily than 
Fic. 2. Lumping of cover sheet area with ribs. Dotted lines 90 , ers i all ion 1j stress 
fall 0.181 of the rib length to either side of rib or midway between ©XP€¢ ted to cause much chordwise bending. sae 
ribs in those cases where the dotted lines would otherwise overlap. The contribution of the sheet to the bending stiffness ol pri 
The bending stiffness of each rib is computed assuming the sheet ‘ 3 eee twee! 
between dotted lines acts with the included rib. of spars, such as the 1 -ading edge spar F, Fig. 3, must be 


torsic 











tlues 9) 
ight, ne 
nsvers 
€ direc 
Oisson's 
3S, both 
in com. 
f delty 
for the 


forma 
d to be 
r sheet 
bs and 
-arried 
V each 
of the 

The 
ONS, is 
otions 
Ctions 
nding 
IS giv- 
of the 


Wise 

with 
bend- 
g the 
ly ef- 


hord- 
sheet 
The 
eness 
r side 
vhere 

rib 
1son- 
their 
toa 
is of 
efer- 
ility 

ap- 
that 
thus 
the 
spar 
r as 
- for 
ily 


1eSS 
t be 








STRUCTURAL 


taken somewhat differently than for the ribs. The 
proper choice of width of flange to be considered as ef- 
jective in a wide flange beam having the flange all on one 
side of the spar web depends on the manner in which 
the moment varies along the beam. As developed in 
reference 9, pages 156 to 161, a reasonable effective 
width to use is 0.181 of the beam length between points 
of zero moment. In the case of ribs, this results in the 
choice discussed earlier and shown schematically in Fig. 
? In the case of a leading-edge spar where substantial 
loads may be applied by the tips of each of the spanwise 
spars, it seems probable that the beam length between 
points of zero moment is more nearly equal to the dis- 
tance between intersections of the leading-edge spar 
with the spanwise spars. This basis for selecting the 
width of sheet acting with a leading-edge spar is pre- 
sented in Fig. 3. It results in considering a minimum of 
sheet as acting with the leading-edge spar and therefore 
has the advantage of considering, again, only a mini- 
mum amount of the sheet that has already been con- 
sidered in connection with the ribs and spanwise spars. 
This method of taking account of a leading-edge spar 
involves the same errors as were noted for ribs. It is 
felt that ordinarily such a spar must be fairly thin for 
aerodynamic reasons and that therefore these approxi- 
mations will not cause serious error in the wing analysis 
isa whole. The errors that result from the probable 
lack of orthogonality between a leading-edge spar and 
the other spars (see, for example, Fig. 3) will tend to be 
small because the stresses in the other spars, at their 
intersection with the leading edge spar, are low. 


Where cutouts are present in the sheet, no considera- 
tion need ordinarily be given to shear lag effects in the 
initial analysis, since it is felt that these will not have 
much effect on the load distribution. In estimating 
maximum stresses, however, particularly in the case of a 
wing root where a spar is clamped while the cover 
sheet is free, it will be necessary to compute the stresses 
on the assumption that, at the wing root, the spar alone 
carries the load that the analysis predicts for the com- 
bination of spar and cover sheet. 


Torque will be considered to be carried by individual 
torque boxes formed from two adjacent ribs, two ad- 
jacent spars, and the thin cover sheets. Each such 
torque box will be considered separately. Its torsional 
stiffness will manifest itself in trying to keep its four 
corners, 7, s, 1, j7, Fig. 4, in a plane. In computing the 
torsional stiffness of these boxes, it is convenient to as- 
sume that the sheet thickness for the spar is infinite. 
This assumption will give little error where (1) adjacent 
torque boxes have nearly equal and opposite shear flows 
in the spar web and (2) the spar web is much thicker 
than the cover sheet. Following the assumption that 
stress along the spars and shear normal to the spars are 
of primary importance, an axis X-X, Fig. 4, midway be- 
tween adjacent spars is taken as the torque axis. The 
torsional stiffness, G/, of the cross section normal to 
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Fic. 3. Lumping of cover sheet area with leading-edge spar. 
Dotted line falls 0.181/;, 0.181, ete., to the side of spar F 


the Y-X axis at its mid-length is then computed as 


representative of the box. 


FLEXIBILITY AND STIFFNESS MATRICES FOR SPARS AND 
RIBS 


The flexibility of the spars and ribs will, in each case, 
first be computed on the basis that one end is clamped, 
as a cantilever beam. This basic computation is then 
modified, as will subsequently be shown, for those 
cases where the end is either simply supported or free. 

Several systematic procedures have been developed 
for computing the bending deflections of tapered canti- 
2 From each of these, the influence 
that is, the deflections at each station due 
can be computed from 


lever beams. 
coefficients 

to unit loads at each station 
the given values of // at each station and the distances 
between stations. The method in reference 12 is 
particularly suitable for use with high-speed electronic 
digital computing machines. The influence coefficients 
in bending for a tapered cantilever beam can therefore 
be considered to be determined by one of these methods. 





.— GJ computed for mid-length 
\ cross-section 


Rib ~< 
Rib 











Fic. 4. Representative torsional stiffness GJ of an individual 
torque box computed for cross section normal to box axis at its 
mid-length on assumption that spar web sheet is infinitely thick. 
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For a particular spar or rib, the influence coefficients 
will be denoted by the matrix [6,;], where 7 and j are 
stations along the spar or rib. Thus, for the tapered 
cantilever beam, the deflection d; at station 7 due to 
RCS Eso Pesy. 5 <3 thas s: « shag AE SUAUIONS.T,.S,.6 J; ss HS 


dy = telly - Cyl + .2. Obey Pc. Baglin (2) 


In the case of a spar with a built-in root, the deflec- 
tions y,; of the spar in the wing are given by the similar 


formula, 


Vi = OipLye + Sighs + ... OjLy + ... binkn (2) 

In matrix form, 
byt = fafn} (2a) 
Lt = [6]—'fy} (2b) 


(clamped root condition) 


In the case of a spar with a simply supported root, 
the deflections y,; of the spar in the wing are, in addition 
to those computed from Eq. (1), those due to the slope 
6 of the spar root. If /; is taken as the distance from 
the spar root to station 7, 


Vi = 10 ca Culy - Bish + see 6,;L; + tee Sinkon (3) 

To eliminate @ from Eq. (3), we make use of the fact 
that the moment is zero at the root 

0 = LL, + IL. + .. bby + .. babe (4) 


/, are the distances, respectively, 
In matrix 


i et a ee) fe 
from the root to stations 7, s, ... 


We 555 Ts 
form, all the equations of the type of Eq. (3) for the 


station deflections can be written 

ly — 10; = [6], Lj (2) 
Inverting this equation gives 

{L} = [6]*'}y — 16} (6) 


{n matrix form, Eq. (4) is written 


Substituting Eq. (6) in Eq. (7) gives 
0 = [/][6]-'}yv — 16} (S) 
Solving Eq. (8) for @ gives 


I 


0 = [5] -'hyvt (9) 
maygy 


Substituting Eq. (9) in Eq. (6) gives 
{7 
{L} = (6 1 BET] | yy (10) 
- Welty 
(simply supported root condition) 
Eq. (10) cannot be inverted to give the flexibility ma- 
trix, since the displacements are not uniquely defined by 
the load but depend on the root slope as well. 
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In the case of a rib or spar having a free root, the de. 
flections y; are, in addition to those given by Eq. (1), 
those due to a displacement y, of the root and a rotation 
6 of the root. 6 may be eliminated as was done for the 
simply supported root, giving [see Eq. (10) | 


oy Pe BAI Tis 4 
i die G May Jeo 


(free root condition ) 





° ) ° ; . 
where the matrices |y — y,{ and }L{ as given in Eq. 
(11) inelude all stations on the rib or spar except the 
root station, ¢/. For a free rib or spar, the sum of the 


forces is zero, or 


—L,-L,—..L1;—... -L (12 


(free root condition ) 


L, = 


When Eqs. (11) are substituted in numerical form into 
Eq. (12), it is seen that the resulting equation, together 
with Eqs. (11), gives all the loads in terms of all the de- 
flections. The flexibility matrix for a free beam can 
not be given, since the displacements are not uniquely 


prescribed by the loads but depend on the slope and dis-’ 


placement of the root. 
For a spar with a root that can be displaced but for 
which the root slope is zero (as would occur when the 


wings are loaded symmetrically ), 
... - =. gener t 
iv — My = [6] L5 


(Ly = [8]-"ty — 5 


(free displacement, zero root slope ) 


(lla 


where }L} and }y — y,} include all stations but the 
root station /. Since, in symmetric motion, the sum of 


the vertical forces is zero, 


(12a 


h=-L-L-..—-L-—.. —-L, 


(free displacement, zero root slope) 


Eqs. (lla) and (12a) give all the loads in terms of all 
the deflections for a spar that can have displacement 


without rotation of the root. 


FLEXIBILITY AND STIFFNESS \ATRIX 
FOR TORSION Box 


The flexibility of the torsion boxes will first be con 
sidered using the plane containing the axes Y-Y and 
)-], Fig. 4, as a reference. Deflections with respect 
to this plane will be denoted by d and loads normal to 
this plane by L’, with suitable subscripts to denote 
panel points’ (superscript 7° denoting the load carried 
by torsion). The rate of twist of the box will be taken, 
see Fig. 4, as 


(L;7); — L,'b ) GJ (13 


where the numerator is the torque carried and the de- 
nominator is the torsional stiffness of the box cross sec- 
tion at the Y-¥ axis, Fig. 4. 
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The deflections are then 


Lim EM 


d; = bx; 
GJ 
L,'b, — L,"b 
d,; = — cb; 7 
GJ 
(14) 
Lh ~ £2," 6, 
d, = —c,b, 
GJ 
Lb, — L,*6 
d, = € b. 


GJ 


The actual deflections yj, yj, v7, ¥s of the box are ob- 
tained by adding the displacements of the reference 
plane to d;, dj, d,, and d;, respectively. Denoting the 


displacement by A and the rotation about the two axes 


Vi lcj(b, + 6.) + o(b, — b;) + c,(b, 
yr feb. — b)) — ob) +b) — €(bi + 6))] + ys[ci(b, + b;) + €)(0i — 0) + cb; + 6 )}) 


} (b, + b,)(b: + b;) (cic; + Cres) + (bi + Os) (0, + 55) (C,€5 + Ci€s) +(b, — b,)(b; — b,)(e:c; + ese )} 


Solving Eqs. (16) and (17) for the loads gives 


~ 


<i 


~s 
& 


\ 


és SS 


where 


"= F(y)[—e(b, + 6.) — (bi + Os) + €s(bi — b,)] 


and }-Y by @ and ¥, respectively, positive when 


they increase the deflection at station 7, we have 
vi = A + 0B, + Wey + Dic ( Li"; — L;"b;)/GI 
yy = A — 0b) + Wey — bye(Li"b, — L,"b,)/GI ; 
ye =X + 0b, — Wor — bc (L,'b; — Lj"b;)/GJ — 
Ys = X — 0b, — Wes + bcAL hh — L,*b)/GI 
From conditions of equilibrium, we have 
Li’ +L,'+L,'+L," =0 
Lb = Lt + Li ~Li, = 68 (16) 


Lia + L;'¢; — Lc — LG = @ 


Eliminating A, 6, and y from Eqs. (15) and solving for 
(L,'b; — L,;'b))/GJ gives 


+ bj)|] + yvil—ci(b, + bs) — c(d) + 0.) + €5(0) — 6,)] + 


‘= F(y)[e(b, + 6.) + (bs — b;)) + c.(b, + 6;)] 


(1S) 


| = Fly)[e(b, — b)) — €(bi + 6.) — es(bi + 5))] 
J = Fly)ledb, + b)) + ¢(b; — b-) + (bi + 5) 


GJ} vile (b, + b,) + c-(b, — bj) + cs(b, + b;)] + yi [—ci(b, + bs) — Cr(b; + bs) + c.(b) — 0,)] + 


V, [ci(d, 


F(y) = 
+ b.)(b; + b)) (cic 


1 (b 


and c = CG = Cj = C, 


When 6 = 6, = 6; = 0, = db, 


i 


The inverse of Eqs. (18), or flexibility matrix, cannot be 
ments; the displacement and rotations of the reference 


INFLUENCE COEFFICIENT MATRIX OF 


COMPOSITE STRUCTURE 


STIFFNESS AND 


The load carried by the composite wing structure 1s 
the sum of the individual loads carried by the spars, 


ribs, and torsion boxes separately. Denoting by P,, 


P., ... P;, ... P, the external loads applied to the wing 
at panel points 7, s, ... j, ... m, respectively, we have 
P= L5+L" +L,’ 
P,=L5+L," +L,’ 
ect ei. oF (19) 
1 


| a ge eee 


— bj) —c(b; + b&b) — Cc 


(db; + b;)] + velco, + 0;) + €)(bi — 67) + (bi + by) ]} 
+ oc.) + (db; + b.)(b, + b))(e¢; + Ces) + (bs — 6;) (bi — b,)(cicy + CsC;) } 


} (d; + b,) (bic, + bcp) + (b, + bs) (bie; + bjc;) 4 


L.7 = Gin — 9 — y+ 
Sh*c =(1Sa) 


-Lj/ = -L,’ = 


Vs) 


given, since the loads do not uniquely define the displace- 


plane are needed in addition to the loads. 


where superscripts S, R, and 7 denote loads carried by 
the spars, ribs, and torsion boxes, respectively. The 
L's in Eqs. (19) can be replaced by their equivalents in 
terms of deflections using, for spars or ribs, Eqs. (2b), 
(10), or (11) and (12), depending on whether the root is 
clamped, simply supported, or free, respectively, and 
using, for torque boxes, Eqs. (18). When this is done, 
Eqs. (19) can be written in matrix form as 

{P} = [A]fy} (20) 
where [A] may be referred to as the composite stiffness 
matrix of the wing. The deformations for a given 
loading condition are given by 


-. 


fy} = [A]-"} Pj (21), 
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The matrix [A]~! is the influence coefficient matrix of 


the wing. 


STRESSES IN WING 


The loads on the individual ribs, spars, or torsion 
boxes can be computed from the deflections given by 
Eq. (21) using Eqs. (2b), (10), (11) and (12), or (18), 
The stresses can then, in 
For in- 


whichever is appropriate. 
turn, be determined by conventional formulas. 
dividual ribs or spars, the loads are used to compute 
bending moments J/ which in turn are used to give stress 
as Mc/I, where c is the distance of the extreme fiber 
from the neutral axis and / is the moment of inertia of 
the section. For torque boxes, the loads are used to 
compute the torque 7 = L,'b; — L,'b;, and from this, 
the shear stress in the cover sheet is taken as 7'/2/A 
where A is the cross-sectional area of the box at the 
)-Y axis, Fig. 4, and ¢ is the sheet thickness. Where 
cutouts and end conditions introduce shear-lag effects, 
the stresses are computed from the loads on the individ- 
ual spars or ribs, taking shear lag into account. 


SOURCES OF ERROR 


Numerous approximations have been made in the 
foregoing analysis in order to obtain a method suitable 
for machine computation which takes account of the 
major bending, twisting, and warping effects that are 
expected in thin, low aspect ratio wings. No account 
is taken of the flexibility of fastenings, although this 
may be significant for heavily loaded, thin wings. 
The effect of Poisson's ratio in producing anticlastic 
curvature has also been neglected. The cover sheet 
has been considered to be thin, although in some appli- 


cations it may be relatively thick. The seriousness of 


AERONAUTICAL 
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these and other sources of error can best be determined 
experimentally. 
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Bending of Thin Elastic Plates of Variable 


Thickness 


* 


Y. C. FUNGt 
California Instilule of Technology 


ABSTRACI 


Southwell’s equations for the bending of thin elastic plates are 
derived from a new point of view and are extended to plates of 
variable thickness and mixed boundary conditions. The bend- 
ing and twisting moments and the transverse shear in the plate 
ire expressed in terms of two “‘stress functions” U and V. The 
boundary conditions are such that, on an unsupported edge, the 
boundary values of U’ and V are known and that, on a clamped 
edge, the first derivatives of U, V are known. Thus, the new 
formulation is particularly adaptable to plates with free edges, for 
Kirchhoff-Love formulation expressed in 
function is inconvenient the 


which the classical 
terms of the deflection 


boundary conditions involve the second- and third-order deriv- 


because 


itives. 

The theory is applied to a square plate with linear thickness 
variation (cross sections in the form of a double wedge) simply 
supported at two diagonal corners and loaded by a pair of equal 
forces at the other two corners. This loading condition is chosen 
for its simplicity in experimental verification that is carried out. 
The stresses and deflection of this plate are found by the relax- 
ation method. The profound influence of the thickness variation 
on the stress distribution is shown by comparing the results with 
the exact solution for a plate of uniform thickness under the same 
loading condition. The theoretical results agree with the ex 
periments within the estimated error of both, indicating that the 
formulation of the bending problems for plates of variable thick 
ness is sound. 

The influence of the thickness variation and the power of the 
method are further illustrated by the results for the bending of 
the square plate by uniform edge moments and also for a full- 


span solid swept wing 


NOTATIONS 


Referring to Fig. 1, the coordinate system, the stress resultants, 


ind the stress moments are defined as follows: 


X,Y, 2 = right-handed rectangular coordinate system 
Vr, My, Hry = resultant bending and twisting moment per 
unit length in plates. The bending mo 
ments are considered positive when they 
put the positive side (outer normal in posi- 
tive c-axis direction) of an element of the 
The 


twisting moment is positive when the in 


plate in tension (positive stresses) 


duced shearing stress rry is positive on the 
positive side of the element. (This con- 
different that used by 


vention is from 
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& 


Timoshenko and is more consistent and 


convenient in the present analysis ) 


Vr, Qy = resultant shearing forces per unit length in a 
plate, normal to the middle surface. Their 
sense is the same as the corresponding sense 
of the shearing stresses, tr2 and ry2, respec- 
tively 

gq = intensity of distributed load per unit area of 
the middle surface of the plate, positive in 
the positive z direction 

Q), Qe = two arbitrary functions defined by Eq. (4) 

x = astress function defining the shear in the plate, 
Eq. (6) 

U,V = stress functions defining the moments in the 
plate, Eq. (8) 

w = vertical displacement of the mid-surface of the 
plate 

D = the flexural rigidity of the plate, = Ef? 
[12(1 — pw?)] 

m = Poisson’s ratio of the material; yp is taken as 
().32 throughout this paper, a value ap 
propriate for aluminum alloys 

E = Young’s modulus of the material; also, in 
Eq. (11), the strain energy 

5, n = elements of length, along and normal to the 
boundary. The positive sense of the outer 
normal m is to that of s as the x-axis is to 
the y-axis 

l,m = the direction cosines of the outer normal to 
the boundary (see Fig. 2) 

On = 1Q, + mQ, = shear per unit length acting on 
the boundary 

MM, Ms. ten bending moments in sections along the bound 


ary, perpendicular to the boundary, and 
twisting moments acting on the boundary, 
respectively, with signs determined by the 
from M,, My, 


tensor transformation law 


and Hy, 


(1) INTRODUCTION 


Shes CLASSICAL KIRCHHOFF-LOVE THEORY of thin 
plates when applied to plates of variable thickness 
results in a fourth-order partial differential equation 
with variable coefficients which is difficult to solve. 
So far, only three cases are known with numerical cer- 
tainty—namely, the symmetrical bending of circular 
plates, the unsymmetrical bending of a circular plate 
with quadratically varying stiffness, and a rectangular 
plate with linear stiffness variation in one direction and 
simply supported on two opposite edges. The sym 
metrical bending of the circular plates has been dis 
cussed by Holzer, Pichler, Olsson, and Conway; the 


unsymmetric bending, by Olsson.'~> The rectangulas 
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plate was also investigated by Olsson, with some addi- 
tions by Reissner.* 7 These solutions are analytic, 
generally involving special functions such as the ex- 
ponential integrals and the confluent hypergeometric 
functions. Earlier mathematical investigations of Birk- 
hoff’ and Garabedian? went deeply into the funda- 
mental formulation of the problem but stopped where 
“the only difficulties involved are apparently in the na- 
ture of disagreeable computation”’ (quote from refer- 
ence 9, page 396). 

For plates with free edges, except those having two 
opposite edges simply supported, few solutions are 
known even for plates of constant thickness. Numert- 
cal solutions have been tried now and then with doubt- 
ful success because of the difficulty of accounting for 
the free-edge boundary conditions, which involve the 
second- and third-order derivatives of the deflection 
function. An efficient procedure, however, was pro 
posed by Southwell, first in 1941" '! and later in 1948." 
Southwell formulated the problem of bending of thin 
plates in two different ways, analogous to the two well- 
known formulations of two-dimensional plane stress 
(or plane strain) problems. One method uses the 
plate deflection as the unknown, and the other uses 
two functions that correspond to the displacements in 
Each type has its special 


the plane stress problems. 
advantages with regard to certain boundary conditions. 
Southwell solved several examples with the relaxation 
In this paper, Southwell’s equations will be 


method. 


SCIENCES JULY, 1933 
derived from a new point of view and extended to plates 
of variable thickness. The functions l’ and 1’ intro. 
duced by Southwell will be identified as two “‘stres 
functions” from which the bending moments and the 
shear in the plate are derived. Practical solutions ar 
obtained by the relaxation method. Since it is a ny 
merical procedure, it can be applied with ease to any 
plate plan form. This is in marked contrast to the 
analytical method, which for mathematical reasons 
suffers from definite restrictions on the geometrical 
shape of the plate. 

It will be shown below that the nonuniformity of the 
plate thickness has profound influences on the stress 


distribution. 


(II) GENERAL FORMULATION 


For a flat plate in static equilibrium, the stress re 
sultants Q,, Q, and the stress moments .1/,, \/,, I, 


(see Fig. 1) are related by the equations 


(OM,/0,) + (Of1,,/O0v) — OQ, = 0 (| 
(OHf,,/0x) + (OM,/dyv) — Q, = 0 (2 
(0Q,/Ox) + (0Q,/Ov) + ¢q = 0 (3 


There are three equations for five unknown fune 
tions M,, M,, H,,, Q,, and Q,. 
these unknowns in terms of two stress functions so that 


It is possible to express 


Eqs. (1), (2), and (3) are identically satisfied. 

When g(x, y) is an analytic function, it is always pos 
sible to find by an integration two functions, 2; and %, 
such that 

g = (0°72), 0x?) + (07Q2/ dv") (4 

In fact, one may put either Q, or 2 to zero, thus clearl 
showing the ease with which Q; or 2, can be obtained 
Here both 2; and Q are retained for the sake of sym 
metry in the formulas that follow.* 


Equation (3) now takes the form 


“(a +) +2 (a +$%) =0 y 
ox oi ox ov % ov) 


This equation is in the form of a divergence of a 
vector field whose rectangular components are Q; + 
(0Q2,/0x) and Q, + (OM, Ov), respectively. Hence, 4 
stream function x can be introduced, so that if 


* In reference 11, Southwell introduced the auxiliary problem 


of determining w,: V?« = qg/D, u = 0 on the boundary, an 
V2w, = u, Ww, = 0 on the boundary, where V? is the Laplacial 
operator. The true deflection w is then put as w = w; + wy», and 


two stress functions relating to w. are introduced_using O(V *z)/0 
and O(V 2z,)/Ov as predetermined quantities. This can be show! 
as equivalent to putting 2; = Q = Q = w, and demanding that 
Q =(Qonthe boundary. In numerical calculations this introduces 
a complication involving the differentiation of a preliminary 
solution 2. The present analysis circumvents this difficulty, and 
hence, the boundary values of l’ and V can be accurately spect 
fied, and the stresses in the plate can be obtained with bette! 


accuracy 
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QO, + (02,/Ox) Ox Ov | (6) etc., or rearranged, 
O 0 /dy) = — xv) f 
2y F (0% /Oy) (Ox, Ox M, = (OV/dy) —Q, 

M, = (OU/Ox) — Q | 


then Eq. (5) is identically satisfied by x. Eqs. (1) and : (S) 
H,, = —(1/2) [((OV/dx) + (OU/dy)] 
X = (1/2) [(0V/dx) — (0U/dy)] 


(2) now become 


The field equations governing LU’ and J’ are of the 


O 
(M, +2) + (UI1,, — x) -9) 
” ” nature of the equations of compatibility. They can 
be derived from the relations among the stress result- 


e 


ants, stress moments, and the deflection function w. 


O O 
(M, Qo) + (11,, + x) = 0 
Ov Ov : ; sie é : 
‘ If the thickness variation is gradual so that its second 


These are again in the form of the divergence of vec- partial derivatives are continuous throughout the plate, 


Therefore thev are identicallv satisfied if the classical assumptions used in the Kirchhoff-Love 


tor fields. 


one introduces two stress functions l’, I’, such that rhe following expressions de- 


theory are applicable. 
rived for plates of constant thickness apply with suffi- 


M, + Q = OV'/Ov, Hy — x = —OV'/Ox cient accuracy to the variable thickness plates also.'” 


-(1 D)M, = (0?w/Ox?) + w (0?w/dv?) = —(1/D) [(OV dy) — Q] | 
—(1/D)M, = (0?w/Oy?) + pw (O°w/dx?) = —(1/D) [(OU' Ox) — Q] (9) 
(1/D) H7,, = —(1 — p) (O°w/Oxvdy) = —(1/2D) [(OV Ox) + (OU dy) ] 


An elimination of w leads to the following equations 
0/1 OU" l1+y40/l1 Ol 0/1 oV 1+yuoO /l oV O /§h Q, 
(2) +2 (53) -2(8- +8) =< 
ox \D ox 2 dv\Doy ox \D dy 2 dy \D ow ox \D D 


Oo /1 ol l+yn0 /1 OV 0/1 OU 1+noO/l ou 0/2) Q. 
( ? ( ) —# ( v = ( — yp 0 
Ov \D oy 2 @x \D ox Ov \D ox 2 Ox \D oy oy \D D 


(10) 


The case in which D = constant was first obtained by Southwell, who also pointed out that Eqs. (10) are anal 
ogous to those governing the stretching of thin plates (plane stress problem) if ’ and V are interpreted as the com 
ponents of displacement. 

An alternative derivation can be obtained by an application of the complementary energy theorem, which states 
that the true stress state is distinguished from all statically correct stress states by the minimization of the comple- 
mentary energy function. This derivation has the advantage of giving the proper natural boundary conditions at 
once. If the boundaries of the plate are either clamped or unsupported (where external tractions are specified), 
the complementary energy function is identical with the elastic strain energy in the plate, 


iat ££ 
E= / / [M2 + M,? — 2u6M,M, + 2011 + pw) H,,*] dx dy (11) 
J J 1 ~+ ot 


his expression may be regarded as the fundamental assumption for the present theory. It replaces Eqs. ({) but 
can be shown to be equivalent to Eqs. (9). 

If \/,, \/,, 17,,, ete., are derived from two arbitrary functions U’ and |’ according to Eqs. (8), the result will repre 
Substituting .1/,, etc., by l’, lV’ according to Eqs. (8) into Eq. (11), one 


(o— ) (> \(>— ) 
—%) — 2u — Q — 2) + 
ov Ov Ox 


| OV. olv\ 7 
=*( + ) | ava (12) 
= Ox Oy : 


Carrying out the first variation of & with respect to 6U and 6V, integrating by parts, and denoting by / and m 
the direction cosines of the outer normal vector to the boundary curve (Fig. 2), one obtains the following expres- 


sent a system of statically correct stresses. 


obtains 


"ff 1 F/ov ; 
E= ( -%) + 
2%1—pw)J J DL\ov 


sion: 
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ae l jf. fav d1 /oVv l+y01/V. dy... 
ee =f J {12 D & 7 2.) ~ ax D E 21) aie dy D = * oy )| me 
oe ae. hr he ada 
ad atl —%) - 9% -9) + d= + vt 
[n BY 2) — am (22a) E81 4 Mar 


The necessary condition that U and V correspond to the actual state of stress is that 5E = 0 for arbitrary 6U and 
6V. Hence, the integrand of the above double integral must vanish, leading to the field Eqs. (10). The vanishing 
of the line integral in Eq. (13) leads to the boundary conditions either 


6U = 6V =0 ‘al 
or 
ou OV l+u aV = 
7 — 1(Q. — wQ,) = O 
(15) 
oe l+yu /oV . ow 
: (25 a sr) a - 1S ™ | — m2, — wh) = | 


The first set of boundary conditions corresponds to the unsupported edge on which edge tractions are specified, 
and the values of U and V are accordingly specified. [See Eqs. (21) and (22) below.] The second set corresponds 
to the clamped edge on which the deflection and slope both vanish. When the entire boundary of the plate has 
only one simply connected clamped edge, the above boundary conditions are sufficient to give a unique solution 
to the problem.* 

The values of U and V on the unsupported boundary are derived as follows: From the coordinate transforma- 
tion relations, 


OU /os = —m(O0U/dx) + (OU /dy) (16) 
Hence, from Eqs. (8), 


OU/dOs = —m(M, + Q) — lay + x) (17) 


i Ox i 02; 02» 
Xx=xot+ | a ds = xo + / (0. + a) +m (0, a a )| ds 
ra "Ss / OQ, O22 
= x0 + i Q, ds + | ( a +m dy ) ds 


* Other boundary conditions such as (1), w and Ow/Omn taking specific values along the boundary, or (2), w = 0 and WM, = 0 
on the boundary, etc., can also be derived from the energy method. In the general cases, however, the complementary en 
ergy function will include a term in addition to Eq. (11). See Sokolnikoff, Mathematical Theory of Elasticity, p. 286; McGraw-Hill 
Book Company, Inc., 1946 


But 


(1S) 


where Q, is the shear per unit length acting on the edge where /,, denotes the moments acting on sections nor 
of the plate, mal to the outer normal n, .\/, denotes the moment 
acting on sections normal to the tangent of the bound- 
QO, = 10, + mQ, (19) ary s, and //,, denotes the twisting couple acting on the 

boundary. 


Furthermore, the tensor transformation of stresses Hence, Eq. (17) can be written as 


leads to oU/ds = —mM,, — 1H, — lx — mQ 


Therefore, 


= 
i 
II 


’?M, + m?M, — 2lmH,, 


M, m?M, + PM, + 2lmH,, (20) ; —_—_ ‘ 
i. =i, = +e ~ won) lia ati -f ee ee ee Ree 
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BENDING OF THIN 


Similarly, one obtains 
V-VW= / (dM, — mH, — mx + 12;) ds (22) 
J? 


The integration constants Ul) and J’) are unimpor- 
tant, since the stresses in the plate are defined by the de- 
rivatives of U’ and V only. If the plate boundary con- 
sists of a single closed curve, it can be shown that Ll’ 
and V defined by Eqs. (21) and (22) are single-valued. 
However, for multiconnected regions, this may not be 
the case. In the following, the discussions will be 
limited to simply connected regions. 

To show that the value of x» can also be arbitrarily 
chosen, note that, from Eqs. (8), xo = constant implies 


that 


(23) 


(OV/Ox) — (OU/Oy) = const. 


If U and V are interpreted as the components of dis 
placement of an elastic body, then the above equation 
means that the body is subjected to a rigid body rota- 
tion. It is clear that the boundary conditions (21) and 
(22) corresponding to x» are consistent with this rigid 
Now a rigid body motion induces no 
stresses—i.e., it does not affect OU’/dx, OV/dy, ete. 
Therefore, 1/7,, J/,, etc., are not affected by the arbi 
trary choice of xo. 

From Eqs. (21) and (22), the values of l/ and V on 
the boundary are completely defined if the tractions 
are known on the boundary. Hence, if the plate 
boundary is completely unrestrained, the problem is 
with given 


body motion. 


formulated into a_ differential 
boundary values rather than boundary derivatives. 
This is an essential simplification over the Kirchhoff 


system 


formulation when unsupported edges are involved. 

For plates of arbitrary plan form and variable thick- 
ness, the orthodox method of integration seems to be 
dificult. The boundary value problem so formulated, 
however, can be solved efficiently by the relaxation 
method. 

Some very simple integrals of the differential system 
seem to be of interest, as given below. 

Consider a plate bent entirely by forces and moments 
applied on the edges and with no distributed lateral 


loading acting. In such cases, g = 0, and one may 
take Q,} = Q. = 0. Consider a solution of Eqs. (10) 
when D = constant, 
U = Mx, V=0 (24) 
This corresponds to a system of forces, 
iM, =-0; Mi= MM e 
(25) 


HH, =x =Q, = Q, = OJ 
Hence, it is an exact solution for a flat plate bent by a 
uniformly distributed bending moment acting on the 
edges y = constant (Fig. 3a). 

This is intuitively obvious. However, it is less evi- 
dent to see why a rectangular plate with variable thick 


ness in the y sections—i.e., D = D(y) (Fig. 3b)—acted 
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on by a uniformly distributed moment on the edges 
y = constant, cannot undergo uniform bending through- 
For if 


M, = Hy = 0, = Q, = 0 


out the plate. 


and 1/7, = a constant, then we must have l’ = Myx 


uv 


and |’ = 0, which, however, does not satisfy Eqs. (10) if 
0(1/D)/oyv ¥ 0 


If, on the other hand, the plate is of variable thickness 
i.e., D = D(x) one 


of the exact solutions is the following: 


in sections parallel to the x-axis 


US const. f D(x) dx, 


which corresponds to the case of bending by moments 
M, distributed along the edges parallel to x-axis, of 


(26) 


magnitude 


M, = 0U/ox = const. D(x) | 
M, = H, =x = 0, = QO, = OJ 


(27) 


Hence, the bending moment is proportional to the local 
flexural rigidity. For example, if the thickness of the 


plate is a linear function of x, 


t=a-+t bx (28) 
then a uniform bending in the plate is possible if 
M, = const.:(a + bx)’ / 
(29) 


M, = H, = 0, = QO, = Of 


(See Fig. 3c.) 

There is no simple generalization of the solution of 
the type (26) for plates with thickness variable in both 
x and y directions. 

(III) BENDING oF A ‘“DOUBLE-WEDGE™’ SQUARE PLATE 
BY UNIFORM EDGE MOMENTS 


In pure bending of plates on which the edge trac 
tions are specified, the boundary values of Ul’ and | 
are known, and the relaxation method can be applied. 
The general principle of the relaxation method can be 
found in the books of Southwell. The particular de 
tails, relating to the examples given below, can be 
found in reference 13. The first step is always the re 
duction of the differential system into dimensionless 
form. Hence, in the following examples, the typical 
lengths and the typical (generalized) forces are always 
one, expressed in terms of the characteristic values. 

Consider first a square plate of constant thickness 
bent by a uniformly distributed moment on two op 
form the boundary 


posite edges. In dimensionless 


conditions are 
M, = 
H,, = Q, = Q, =q = O0onx = 


= = |] } 
+ | \ 


lony = +1, M, = Oonx 


(30) 
Liy= 
Evidently an exact solution is given by Eqs. (24) and 


(25), provided that J = 1. 
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yfo form, .\/, is everywhere zero and J\/, is constant (= |). 
\w3wwarawwrwrrrree For a double-wedge plate, however, ./, is no longer 
ly Ft to zero, and M, becomes variable. While V/, differs only 
PP eh Reais 2 slightly from unity, 1/, attains values as high as 0,12 
Tj BAATA RE SBA : z J : , : : ; 
7 (1.e., approximately 12 per cent of .1/,) in certain regions 
¢ of the plate. 
~ fo Unfortunately, an exact solution to the above prob- 
y; prob 
& i x lem is not known. Therefore, it is impossible to esti- 
mate the accuracy of the relaxation procedure with re- 
gard to the requirements on the fineness of the nets, 
On the other hand, experimental foundation is also 
€ ££ ££EELE lacking in the theory of plates of variable thickness, 
Although the fundamental assumptions made in de 
FIGURE 4 riving the thin-plate equations have been extensively 





In contrast to this simple result, let us consider a 
square plate of variable thickness with the same bound- 
ary moments as given by Eqs. (30). Let the thickness 


distribution be defined by the equation 
t = %& [1 — (ly}/2)] (31) 
so that the flexural rigidity variation is 
D = Dy [1 — (jy|/2)}* (32) 
where 
Do = Eto?/12(1 — yw?) (3:3) 


The cross section of the plate in the y direction (v = 
constant) is therefore a ‘‘double wedge”’ (Fig. 4). 

The results are given in Fig. 5. The values of LU’, I’ 
and .V/,, \/, are recorded at each nodal point and are 
arranged as shown in the legend. Contour lines of 
equal magnitude of 17, and M, are also shown. 

The effect of the variation of the plate thickness is 


clearly shown by Fig. 5. If the plate thickness is uni- 


tested for flat plates, few experiments have been carried 
out on plates of variable thickness. This lack of ex- 
perimental verification is probably due to the fact that 
there are but few analytic solutions known with nu- 
merical certainty, and, hence, there has been little to 
check against. It is therefore thought worth while to 
work out an example that can be tested easily, thereby 
verifying the fundamental assumptions. This leads to 
the example in the following section. 


(IV) BENDING OF A “‘DOUBLE-WEDGE” SQUARE PLATE 
BY LOADS AT THE CORNERS 


Consider the same double-wedge square plate with 
thickness distribution given by Eq. (31). Let this 
plate be acted upon by two pairs of concentrated 
forces in opposite directions at its four corners. Ac 
cording to the classical theory of thin plates, this load 
ing is statically equivalent to a system of uniformly 
distributed twisting moments //,, acting along the edges. 
Let the origin of the (x, y)-axes be taken at the center of 
the plate and the semichord length be taken as the char 
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acteristic length. One may consider the plate as 
loaded by 

H,, = 1, M, = Q, = 0 (34) 
on all edges. Hence, 2; = 2, = 0, and the boundary 


values of the functions U and V are easily seen to be 
those shown in Fig. 7. 
It can be shown that if 


M, = M, = Hy = 0 


on the boundary but Q, are concentrated at the corners 
so that 

S Q,ds = +2 
if the integral is extended across a corner, U and V are 
reducible to the same boundary values given in Fig. 7, 
except possibly differing by an integration constant 
that has no effect on the stress distribution. 

Because of the symmetry and antisymmetry of the 
boundary values, only one quadrant of the plate needs 
to be recorded. The results are given in Figs. 8-11. 
A net with 81 points (including the centerlines) in a 
quadrant is chosen (189 points for the full plate). The 
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FIG.8 — BENDING AND TWISTING MOMENTS IN DOUBLE WEDGE SQUARE PLATE 7 
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moments M/,, M,, H,, are tabulated in Fig. 8, in which 
the arrangement at each nodal point of the net is the 
following: 


M, 
M, 


H,, 


It is clearly seen that the normal component of the 
bending moment, /,, vanishes everywhere along the 
boundary. //,, does not remain unity on the boundary, 
but it is easily verified that the reaction Q, + (O/7,,/Os) 
Hence, the static condi- 
The values of 


does vanish on the boundary. 
tions on the free edges are satisfied. 
HT,, at the corners of the plate are unity. 
to the classical theory of thin plates, this corresponds to 


According 


a concentrated force acting at the corner with a mag- 
nitude 


F = 2H, = +2 


The contour lines of 7, and M, are shown in Fig. 9. 

From the data in Fig. 8, the principal flexural mo- 
ments 1, MM, and the angle of inclination of the prin- 
cipal axes, a, can be computed according to the formu- 


las 
M,+ mM, (M, — M,)? 
M, += La 
2 4 (35) 
tan 2a = —2H,,/(M, — M,) 
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The resuits are given in Fig. 10. The crosses at each 
nodal point of the net show approximately both the 
magnitudes and the inclinations of the principal mo- 
ments. 


The outer fiber stresses can be easily obtained from 
6M,/? | 


o, = 6M,/?’, oy, = 


= 6H,,/t2, 


(36) 


Try “23 6M. 2 t2\ 


where ¢ is the local thickness of the plate. Using these 
equations, the maximum principal stresses can be com- 
The contour lines of the maximum principal 


The numerical values 


puted. 
stresses are plotted in Fig. 11. 
attached to these curves are the values of 


to” C= 6M to? ‘t? 


where o and MW are the maximum principal stress and 
moment, respectively. 

In order to compare the theoretical results with the 
strain-gage measurements, the outer fiber strains are 


required. These are obtainable from the following re- 
lations: 
(Et?/6) «- = M, — uM,) 
» (37) 
(Et?/6) «, = M, — uM,) 


All the numerical values relating to Figs. 8-11 are 
reduced quantities. In returning to physical dimen- 
sions, if the length of the edges of the square plate is 2L 
in. and the force applied at each corner is 2P Ibs., then 
H,, = P (tb. in. per in.). The bending moments 1,, 
M,, the twisting moment //,,, and the principal mo- 
ments J/;, 7, in the plate (Ib. in. per in.) are equal to 
the values given in Figs. 8, 9, and 10 multiplied by P. 
The stresses o in pounds per square inch are given by 
the numbers in Fig. 11 multiplied by P/ty*, where fo is the 
thickness of the plate at the centerline y = 0, in inches. 

The contrast between the behavior of a flat plate and 
plate of variable thickness is remarkable. If the thick- 
ness is constant, the exact solution to the above prob- 


lem is 
Hay = I / 
} (38) 
M, = M, = Q, = Q, = 0) 
throughout the plate. Hence, for a flat plate, instead of 
Fig. 8, one would have M, = M, = Oand H,, = 1; in- 


stead of Fig. 9, one would have a complete square of 
zeros; instead of Fig. 10, one would have all crosses 
with unit length at 45 deg. inclination everywhere; and 
instead of Fig. 11, one would have uniform magnitude 
everywhere. The differences are so remarkable, and 
at first sight so surprising, that an experimental veri- 
fication necessary. Thus, Fig. 9 shows that 


M,, which in a flat plate would be zero, can become as 


seems 


high as SO per cent of the principal moment for a wedge- 
shaped plate. Fig. 10 shows that, in certain regions 
of the plate of variable thickness, the principal moment 
can become almost twice its value in a corresponding 
flat plate. The angle of inclination of the principal 
axes or of the lines of principal curvature can rotate 
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as much as 15 deg. All these remarkable differences 
occur in a plate with linear thickness variation and ce, san b 
with the thickness at the middle line only twice that P 


at the edge. 

An experiment was therefore carried out to check 
the theory. A 10-in.sq. plate of mild steel was made 
to the dimensions of Fig. 12. Fifteen Baldwin-South- 
wark A-7 electric wire resistance strain gages were put 
on the plate, plus two standard AR-1 strain rosettes 
manufactured by the same corporation. The gage 
locations, label numbers, and the gage constants were 
indicated in Fig. 12. Regular Duco cement was used to 
mount the gages upon the plate. The electrical circuit 
was the same as that described in reference 14. The 
details of the measuring equipment, test procedure, 
and the data reduction methods were also the same. 
Hence, reference should be made to that paper for 
detail discussions on the experimental technique. The 
plate was supported at two opposite corners on two 
small spherical pivots on the table. The third corner 
was held down by another simple support pointing 
downward and rigidly attached to the table. The 
loading was then applied at the fourth corner by pre- 
measured lead weights. Since the loading system was 
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FIGURE II 
THE MAXIMUM PRINCIPAL OUTER FIBER STRESS 
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FIGURE 12 
DOUBLE WEDGE TEST PLATE 
LEGEND 
0.225" STRAIN GAGES: 
oO 50.10" “I. P A-7 
1 7 os =o =o . . ies GAGE DIRECTION 
1.25" fe 1.28" 1.25" 1.25" GAGE NUMBER 
a+ GAGE FACTOR =1.97 
he 1.25" 125" 
oO 
78 { : - AR-1 ROSETTE 
168 -_ 0.250" = GAGE NUMBER 
54 GAGE DIRECTION 
6 GAGE FACTOR: 2.06 
>|} 0.10" 
one 16 B, 218, ETC. - GAGES 
=e ON BOTTOM SIDE OF 
PLATE IN SAME POSITION 
AND DIRECTION 
& e | 10.00" 
——— 260" b~— 2.50" 
ow 
‘ Fi 
la 10.00" ~ 
0 -0.30 ~0.10 -0.02 fe) statically determinate, any support deformations, pro- 
O 0.10 0.03 0.0l fa) vided they were sufficiently small, were of no concern 
to the results. By adding the lead weights one by one 
successively, the linearity of the measuring system 
1.53 009 0.53 0.23 O could be checked. In the experiment, a maximum load 
-4.79|-4.33 -236 -080|-070 -04! O of approximately 200 Ibs. was applied. The linearity 
of the potentiometer readings against the applied loads 
was satisfactory for most of the gages, but gages Nos. 
1.31 102 0.6! 028 O 3, 4, and 5 indicated nonlinearity. Consequently, 
-41! |-3.62 -2.38+193 -I.1| +097 O50 O they were replaced, and another series of tests was 
performed. However, it was not able to get linear 
readings* (i.e., the difference in the potentiometer 
0.90 0.68 0.41 0.19 O readings did not vary linearly with the difference in 
-2.61 |-2.42 -1.60 06! | 0.63 0.35 Oo loading). Hence, the data from these gages were not 
presented. For other gages, the data were fitted with 
the least-squares linear regression lines, and the changes 
fe) O O O o| in millivolts of the potentiometer readings at a load of 
O re) O O O 125.7 Ibs. at each corner were determined. The strain 
was obtained from the formula 
FIG. 13— OUTER- FIBER STRAIN 
COMPARISON BETWEEN THEORY AND 4 (millivolt reading) 2g 
a e= SE Cnn aE (39) 
EXPERIMENT FOR THE GORNER-LOADED (gage factor) - (battery reading) 
DOUBLE WEDGE PLATE : 
THEORY EXP * A recent investigation shows that this nonlinearity is caused 
€.-104 € 104 by the ‘‘membrane stresses’”’ (stresses in the middle plane of the 
x 3 x r plate), which are neglected in the above formulation. This phe- 
€y 10 €y!0 nomenon of a “boundary layer” caused by the nonlinear terms in 


the plate equation is to be published in another paper. 
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where € is expressed in inches per inch times 10~*. The 
results of these tests are shown against the theoretical 
values in Fig. 13. 

Since the gages Nos. 1, 2, 7, 9, 19B, 20B, and 21B 
were located at the same corresponding point in dif- 
ferent quadrants and on opposite sides of the plate 
the strain readings at these points should be equal 
The scattering of the data, therefore, 
For these 


numerically. 
gave some idea of the accuracy of the test. 
seven gages, the standard deviation of the gage read- 
ings was 0.0073. The mean reading of these gages was 
0.0973. Hence, the experimental 95 per cent confi 
dence limits of the strain-gage readings ¢€, at this loca- 
tion (the midpoint of each quadrant) were 0.0973 + 
0.0060. This is to be compared with the theoretical 
value 0.111 in Fig. 13. The theoretical value is com- 
puted by assuming the surfaces of the plate to be per 
fectly flat. The effect of the finite mesh size in the re- 
laxation process is uncertain. Moreover, the finite ex- 
tent of the A-7 strain gages ('/, by 
average value of e, over the area occupied by these 
uncertainties shows 


3/16 in.) gives an 


gages. Consideration of these 
that the difference is really insignificant. 

On percentage basis, the possible error appeared to 
be higher than that claimed in reference 14. If a 
single measurement is made at one point, the strain- 
gage reading may have an uncertainty of +15 per cent 
for 5 per cent significance level. (The tests reported 
here were made by different personnel from those re 
ported in reference 14.) 

Although a tendency of the theory to overestimate 
slightly the outer fiber strain seems evident in Fig. 13, 
all of the experimental and theoretical values of strain 
Hence the difference be 
There 


agreed within 15 per cent. 
tween theory and experiment was insignificant. 
fore there is no reason to doubt the correctness of the 
theory on the basis of these tests. 

From this example it appears that better technique 
in the experimentation is needed in order to check the 
accuracy of the fundamental assumptions made in the 
for variable-thickness 


derivation of the equations 


plates. 
(V) A So_tp Swept WING 


To avoid undue complications, let us consider the 
A wing of the 
sweep angle. 


following relatively simple example. 

plan form shown in Fig. 14 has a 45 
The cross-sectional shape is flat from the leading edge 
to the mid-chord line and then straight taper to half its 
thickness at the trailing edge. The thickness of the 
center panel is the same as that of the outer panel. 
There is no taper in the spanwise direction. The wing 
is attached to the fuselage by four simple pin joints 
symmetrically located at Ri, Re, etc., as shown in the 
figure. Since the aerodynamic loading can be divided 
into a symmetrical and an antisymmetrical part, it is 
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FIG. 14 


DIMENSIONS OF A SWEPT WING 





clear that the reaction on these joints due to the aero 
dynamic loading can be determined statically. (The 
redundant stresses induced by the indeterminacy of the 
joints, as by misalignment, etc., are of no concern in our 
analysis.) Let us consider a case of uniform symmet 
rical wing loading—.e., 


q = const. = q (40) 


throughout the wing surface (outboard of the wing 
root EB, Fig. 14). The assumed 
unloaded except by the 


center panel is 


fuselage pin-joint reac- 
tions. 

The simplicity of the external loading, Eq. (40), is in 
significant in reducing the labor required for the solu 
tion. If the exact aerodynamic loading qg (x, y) is 
known, the computation can be carried out with a 
small amount of additional labor. On the other hand, 
the simple geometry of the wing plan form is an im- 
portant simplification. Otherwise, all the boundary 
of the wing would not fall on the nodal points of the 
relaxation nets, and the “‘irregular stars’’ will introduce 
a serious question of the computational accuracy. 
The treatment of the “irregular stars’’ 
cussed at length by Southwell.'® Since the purpose of 
the present analysis is to show the effect of sweep and 
of variable thickness, an example without irregular stars 


has been dis- 


seems desirable. 

Since g = qo is symmetrical with respect to the fuse- 
lage centerline, 2, W/,, /,, and V must be symmetrical, 
and U’ must be antisymmetrical, with respect to the 
same centerline. This shortens the work by half. 
Only the half-wing is recorded in the subsequent 
The half is determined by sym- 


solution. other 


metry. 
Let the semichord length be taken as the character- 
istic length, and let 


gq = | (positive upward) (41) 


over the entire cantilever wing outboard of the root 
section EB. Then, it can be easily shown that the 


reactions at the pin joints are 
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FIG. 15 
xy IN THE SWEPT WING 


BENDING MOMENT IN THE WING 
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‘ Hence, the boundary values of U,and V can be eval- 
R,; = 1, upward, so J Q,ds = 1 a 3 ™ a 
Ry uated according to Eqs. (21) and (22). 
(42) rhe problem is again solved by the relaxation method. 


Ke = 


(1 
2 


0, 


JI 
x — $y), 


~~ 
~ 


we 


~*~ 
~ 
_ 
— 


—11, downward, so fo. ds 
Re 


For the loading specified above, one may take 
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——— 


can 


about 120 man-hours. From the accepted values of 

—" a /y U and V, the stresses can be determined. The results 
r are given in Figs. 15 and 16. The contour lines of the 

“ (43) maximum principal moments and stresses are also 

forx < Z a plotted in these figures. It is interesting to observe 


A solution for a rough net with 19 nodes in a half-wing 
next re- 
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be carried out in about 12 hours. 


finement, with a net of 86 nodes in the half-wing, needs 


the stress concentration at the wing root of the trail- 


ing edge. 
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In returning to the physical units, if the chord length 
(measured perpendicular to the leading edge) is 2Z in. 
and the wing loading is g lbs. per sq.in., then the mo- 


M,, Hy 


multiplying the numbers in Fig. 


ments MM (in.lbs. per in.) are obtainable by 
15 by gL’, and the 
stresses in the plates (Ibs. per sq.in.) are obtainable 
by multiplying the numbers in Fig. 16 by gL*/to?, where 
ty is the thickness of the wing at the mid-chord point in 
inches. 


(VI) CONCLUSIONS 


Southwell’s formulation of the problem of flexure of 
thin flat plates can be extended to plates of variable 
thicknesses. An application of the relaxation method 
then renders a general treatment of such plates possible. 
In this paper several examples are worked out. These 
examples show that the effect of variable thickness on 
the stress distribution is large. Hence, in engineering 


designs, if one idealizes a variable-thickness plate into a 
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flat plate, great care must be exercised, even if the 
variation in thickness is small, lest serious error result. 

The detail techniques of the relaxation procedure are 
not presented in this paper. A few useful remarks gath- 
ered from the experience of actually working out 
the examples can be found in reference 13. A method 
of checking the boundary values of U, V on the free 
edges is also given there. Since in the present method 
the external loadings are expressed in terms of the 
boundary values of U and V, any mistake in the evalu- 
ation of the boundary integrals (21) and (22) will cor- 
respond to a wrong specification of the external load- 
ing. It is therefore extremely important to be assured 
of the boundary values. 

The evaluation of the accuracy of the relaxation 
method and the determination of a criterion to decide 
the fineness of the relaxation net needed to achieve a 
specified accuracy are impossible at the moment, be- 
cause no exact solution exists for these illustrative ex- 
amples. Comparison could be made with experimental 
results, although this would involve not only the ap- 
proximations introduced by the relaxation process in 
solving the differential equations but also the approxi- 
mations implied by the underlying assumptions made 
in the derivation of these differential equations. Care- 
ful experimentation and accurate theoretical numerical 
solutions should be carried out in the future to justify 
the fundamental formulation of the variable-thickness 
plates problem. 
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Equations for the Loading on Triangular 
Wings Having Supersonic Leading and 
Trailing Edges Due to Various Basic 
Twist Distributions' 


JULIAN H. KAINER? 


Republic Aviation Corporation 


(1) SUMMARY 


The linearized supersonic flow theory has been applied to de- 
termine the potential, pressure and span-load distributions due to 
various basic twist distributions on triangular wings with lead 
ing and trailing edges ahead of the Mach cones. 

Generalized equations in closed form are presented for the po 
tential, pressure and span-load distributions due to constant, 
linear, parabolic, and cubic twist distributions. In addition, de- 
sign charts for the rapid evaluation of the pressure and span 
load distributions are presented. 

The results are directly applicable to rigid wings with built-in 
twist and may be applied to the aeroelastic problem by means of 
superposition of the basic solutions. In addition, a means of 
obtaining pressure distribution on wings in the presence of bodies 


is suggested 


(Il) SyMBOLS 


Free-Stream Conditions: 


V = velocity 

M = Mach Number 

8 = VM? -1 

m = Mach angle, aresin (1/17) 
p = mass density of air 

q = dynamic pressure |('/2)p V?| 


Wing Geometry: 


b/2 = half-span 

( = root chord 

Cav = average chord 

€ = half apex angle, deg. 

6 = angle between trailing edge and axis, deg 

m = Ptane 

m = Btan 6 

Ss = wing area 

X, 9, 2 = Cartesian coordinates of system of axes with 
origin at leading edge of root chord 

&,”, ¢ = integration variables 

a(t, n) = angle of attack in the stream direction, rad 


Analysis Parameters: 


o = velocity potential 

¢’ = elementary solution to linearized equation 

ox = horizontal perturbation velocity 

dsr, Gyy, Gee = second partial derivatives of the velocity po 


tential with respect to x, y, 2, respectively 
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c(é, n) source-strength per unit area 

w = upwash 

f(g, n) = 1/V (x — &)? — By — n)? — Bs — £)? 

Y = an arbitrary characteristic dimension 

n = nondimensional span coordinate 

a = area of integration 

t = By/x 

A, = cos! [(1 — mt)/(m — 2)] 

Os = cos! [(1 + mt)/(m + 2)] 

x = x/¥ 

ad = x/co evaluated at the trailing edge 

Ap/q = lifting pressure coefficient 

(Ap/q)* = effective lifting pressure coefficient per rad 

C, = section normal force coefficient per rad. 

te)" = effective section normal force coefficient per 
rad. 

a = coefficients of power series expansion for the 


angle of attack 


(III) INTRODUCTION 


pen PRESSURE DISTRIBUTION and the section nor- 
mal-force distribution on twisted wings requires a 
knowledge of the distribution of twist. For rigid wings 
wherein the twist is built-in, the problem may be solved 
directly, since the strength of the source function is 
known everywhere. In aeroelastic problems, how- 
ever, the distribution of twist upon which the pressure 
and span-load distributions depend is unknown. A 
powerful method! is described to obtain the twist dis- 
tribution and, hence, the pressure and span-load dis- 
tributions for equilibrium without recourse to exhaus- 
tive iterative procedures.” * This method depends 
upon the knowledge of the loading due to various known 
basic twist distributions. 

The present report, therefore, has a twofold purpose: 
firstly, to present the basic span-load distributions 
from which solutions may be obtained using the method 
referred to above! for the aeroelastic problem, and, 
secondly, to present the means of obtaining the po- 
tential, pressure, and span-load distributions due to 
known basic twist distributions for rigid triangular 
wings having supersonic leading and trailing edges. 

The lift, moment, and center-of-pressure for the pres- 
ent physical problem will be reported on at a later 


date 
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Details and limits of integration. 


(IV) ANALYSIS 


(A) General Background 


The well-known linearized equation of motion in 


three dimensions is 


B*drr — Gy — O22 = O (1) 
a solution of which is‘ 
o’ = c(é, n)f(é, n) (2) 
where 
f(é 0) = 


(1/V (x — £2? — By — n)? — Bis — ©)*)T (3) 
and c(£, 7) is the source strength per unit area, shown? 
to be 


= (V/mr)alé, n) (4) 


. 


c(é, n) = w 
Within the limitations of the linear theory, it is per- 


missible to separate the effects of twist and camber 


a(é, n) = a(n) + a(é) (5) 
and we are concerned only with the effects of a(n) herein. 
The twist distribution, a(n), is expanded in a power 
series about the root chord 


a(n) = Ya (n/¥)° (6) 


i=0 


)2 is not considered, inasmuch as 


+ Hereafter the term 6%(s — ¢ 
the solution must be satisfied in the plane of the wing. 


“JULY, 196s 


where the absolute value signs imply symmetry about 
the x-axis. 


(B) Velocity Potential 

The velocity potential at a point P(x, y), Fig. 1, js 
the integrated effect of the sources distributed within 
the fore-cone of P; hence, from Eqs. (2) and (4), it 
follows that 


V t 
@ = R.P. | LS a(n) f(é, n)dé in|! (7 
© JS: 


where R.P. means the real part of the expression. 

Upon substituting the twist distribution, Eq. (6), 
into the equation for the potential, Eq. (7), it follows 
that 


bd = Ao(/a)o + A1(o/a), + A2(G/a)o +... (Sa) 


or 


od = 2. a,(o/a); (Sb) 


7=0 


because the integrand resulting therefrom is absolutely 
convergent—4.e., 


iD ai\(n/y)‘|f(E ) < YL iad(n/y)|fE 0} (9) 
7=0 7=0 


The potential for any basic twist distribution is 


Vy -7 
(o/a); = J J In'| f(é, n)dé dn (10) 
| s! 


integration of which leads to results of the form 


t.,% * 
(7) (22) (*)°= aman 
J git? Ney 


The right-hand member of Eq. (11) is presented in the 
Appendix for values of ¢ = 0, 1, 2, and 3. 


(C) Pressure Distribution 


The lifting pressure distribution is defined as 


Ap/q = (4/V) ¢; (12a) 
| o *) 
= he (12b) 
V p> Ov (° i 
= > a,(Ap/qa); (12¢) 
7=0 


From Eq. (11) the horizontal perturbation velocity 


parameter is obtained as 


\ (8\' #)° ) | 
od A(t) —t A,(t (13a) 
(=) (®) (¢ (1 + Ai > | i(t) | 1 


= B,(t) (13b) 
which is presented in the Appendix for values of 7 = 0, 
1,2, and 3. 


t In all the subsequent formulas involving potential, pressure, 
and span-loading, the real part will be implied but the symbol 
R.P. will be deleted 
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Conical Kay, t= 44x 


Fic. 2. Pressure coefficient parameter, (Ap/ga)o*, against conical ray, t, for several families of 8 tan ¢ due to constant spanwise twist 
I q 8 . . ° 
distribution 


D) Span-Load Distribution functions represents a check on the mathematical der- 
ivations, they are presented in the Appendix as A,(t) 


The span-load distribution, or section normal-force 


coefficient, ¢,¢ is defined as and B,(t), respectively. 


f E Ap F 1) 
CzcC = ax ( 2 
ZL. q 10 


Substituting Eq. (12b) in Eq. (14) results in 


| “TI oO 
= — | » # a, (°) dx (15) 
| T11 er ox \a 


and upon consideration of the absolute convergence of 





S) 


the integrand 


17) 

= 
oodr7opoz 
Oo 


BYX 


OC = > a;(@/ a) (16) 


$ =( 


where (¢/a); is evaluated at the trailing edge, since the 
potential vanishes at the leading edge. 

A convenient span-loading parameter may be con- 
structed with the aid of Eq. (11) in which the char- 


acteristic length, y, is taken as the semispan; hence, 


a mm, \* /1\'!+t# fe,\ * c _ 
a = A,(t) (17) 
S Myo — Mm x Gis Con 


where / and x’ are evaluated at the trailing edge. It is 
pointed out that A ,(/) is presented in the Appendix for , 


¢= 0, 1,2. and 3. LO 


Conical Ray ,t 





10 10 10 
E) Unswept Wings j Ci 
nswept Wings Lifting Pressure Coefficient 
rhe limiting cases for the functions AK!) and B,(t) Parameter, (4/X)(4PAY) 
are obtained when 7 = ©, which implies unswept 1 


wings (infinite Mach Number is bevond the scope ot the HIG » Lilting pre ssure reels Incient parame ter, (9 x (Ap qa < 
: : igainst conical rav, ¢t, for several families of 8 tan e due to linear 


linear theory). Inasmuch as the limiting case for these symmetric twist distribution 
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Fic. 4. Lifting pressure coefficient parameter, (8/X)? (Ap/qa)2*, against conical ray, ¢t, for several families of B tan « due to parabolic 
symmetric twist distribution. 


(V) DiIscuSSION 


(A) Load Distribution 


The pressure distributions and the span-load dis- 
tributions resulting from the present analysis are pre- 
sented in Figs. 2 to 5 and 6 to 9, respectively, for several 
families of 8 tan e and with regard to the span-load distri- 
butions for all families of 6 tan 6 satisfying the triangular 
boundary condition (supersonic leading and _ trailing 
edges). The span-load distribution parameter is a 
function of the trailing-edge sweep; however, it is inter- 
esting to note that the results of Figs. 6 to 9, which were 
calculated independently of the trailing-edge sweep, 
are in reality valid for all 8 tan 6! This artifice was 
obtained through the use of the section normal-force 
coefficient parameter plotted against the conical param- 
eter = By/x from the root chord (f = 0) to the leading 
edge (f = m) instead of plotting the section normal- 
force coefficient against semispan, y’. Since the fol- 
lowing equations relate x’ and y’ in terms of mp and ¢, 


it is apparent that the results of Figs. 6 to 9 may be 
used to obtain the span-load coefficient against y’ for 
all supersonic values of 8 tan 6! 


x’ = mo/(mo — 2) (18) 


y’ = (m — mo)t/m(my — 6) (19) 


This method of calculating the span-loading param- 
eter has further merit insofar as the equations for the 
pressure distribution have similar functions of f 
namely @,, 6, and V1 — f*—although in the former 
case / is based on the value at the trailing edge, whereas 
the latter is more general. However, this does not 
matter. 


(B) Application 


In applying these results to the aeroelastic problem,’ 
it should be pointed out to the reader that the method 
involved therein’ is not iterative. Although each basic 
twist distribution gives rise to an elastic deflection that 
may be represented by a power series, it is the cceffi- 
cients of these power series which reduce the problem 
to one of a certain number of unknown constants with 
a compatible number of simultaneous equations. 

Furthermore, if the internal structure is of such na- 
ture that the elastic deflections are irregular, the 
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— 
© 


method! and solutions given herein may be well adapted (VI) ConcLusIons 


to matrix curve-fitting techniques.® oe s 

- (A) The results of the present paper have been pre- 
2) Whine Badly dnote sented in such a form where they are applicable to (1) 
: ; the direct problem of determining the pressure distribu- 


Beskin’ showed that the upwash caused by a body of — tion and load distribution on a triangular wing given 


revolution at a finite angle of attack in a supersonic 
stream is maximum at the body and decreases asymp- 
totically with distance from the body. Essentially, 
this merely distorts the air stream in some man- 
ner. 

If a rigid wing is placed in such a distorted stream, a 
conceivable analogy is that of a twisted or twisted and 
cambered wing in a uniform stream. Since the effects 
of twist and camber are separable, the results presented 
herein might lend themselves to empirical improve- 
ment of the present means of predicting span-loadings 
on wing-body combinations. Experimental verifica- 


tion may be justified. 


the twist distribution, and (2) the inverse aeroelastic 
problem wherein the twist distribution is unknown. 
In the latter case, the twist distribution is determined 
by a method described! from which the final loading is 
determined with no recourse to iteration. 

(B) The span-loading is presented in a concise form 
independent of the trailing-edge sweep but is valid for 
all supersonic trailing-edge sweeps within the scope of 
this report. 

(C) It is recommended that the value of using these 
results to obtain a higher order approximation to the 
problem of wing-body interference be experimentally 


investigated. 











(VII) APPENDIX 


The conical functions that prescribe the velocity potential, the pressure distribution, and the span-load dis- 
tribution for 0, 1, 2, and 3 degrees of twist for triangular wings having supersonic leading edges are presented in 
this section. For the sake of completeness and for checking purposes, the corresponding results for the limiting 
case 8 tan e = © (unswept wing or infinite Mach Number) are also presented. 


(A) General Triangular Wings 


A,(t) = (1/V m? — 1) [(m — DO + (m + 1A] (20) 


A,(t) = son? Vm? —-1V1—-2 + (?(m? — 1)” arecosh (1/t) — 


l l 
= [m — (2m? — 1)t] (m — tI) — = [mz + (2m? — 1)t]}(m + 10st (21) 
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] l 
A.(t) = —e zr. J [(m? + 2) (m — t)? — 6m(m? — 1)t(1 — mt)] (m — 16 += X 
3(m? — 1)” \2 2 


[(m? + 2) (m + 1)? + 6m(m? — 1)t(1 + mt)| (m + 1 — m Vm? — 1 (3m? + (5m? — 2)t?] vi ~ re (22) 


3 V m? — 1 


re (1 — 1llm? — 10m‘*)t? V1 — {?+4+ 
2m? 


n4 1 
Ax(t) = | om }~ (4m? + 11) Vm? — 1(1 ~ oy 


(m2 — 1)” YQ 
, l 3 

: (m? — 1)’*t4 arccosh (1/t) — 3 F (2 + 3m?) — m(4 + m*)t + “(1 + 4m?)t? — m(3 + 2m?) — 
m 2 


l 3 
(2 — 7m? + 8m* — 8m'*)t* | 0, — 3 (2 + 3m?) + m(4 + m?*)t + : (1 + 4m?)t? + m(3 + 2m?)t, — 
4 2 


4m*4 
l se ; ) 
- (2 — 7m? + 8Sm* — Sm®)t* | A (23) 
4m?4 f 
By(t) — ee. 4 
0 = = 1 2 (24) 
Vm? — | 
m? 
B,(t) = —ioney Wert ~ 191-97 -— C0 + mt) 0, — (1 + mt)@]} (25) 
(m2? — 1)” 
m3 / 
B,(t) = 7 1” t [(m — t)? + 2(1 — mt)?]6, + [(m + 1)? + 2(1 + mt)?|& — 6 Vm? —1V1 — #F (26) 
zim? — ex 
m4 / ee / : 
B;(t) = {2 Vm? — 1 [4m? + 11 + (4 + 11m*)t??] V1 — ? - 3(1 — mt) X 


6(m? — 1)” 
[2(1 — mt)? +3(m — t)?]0, — 3(1 + mt) [211 + mt)? + 3(m + £)?]@} (27) 


(B) Unswept Wings 


A(t) = (28) 
A(t) = [V1 — t? + 2t arcsin (t) + ¢? arccosh (1/f)] (29) 
A(t) = (w/6) (1 + 68?) (30) 
= ae /- 

’ A;(t) = - Fe + 13f) V1 — # + # arcecosh (1/t) + 22(1 + 22°) aresin o | (31) 
By(t) = 7 (32) 
B(t) = 2 [V1 — t? + ¢ arcsin (t)] (33) 
B,(t) = (m/2) (1 + 2) (34) 
’ r . | . ; 
B;(t) = : | + 11?) V1—# + 3t (3 + 2°) arcsin | (35) 
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Review of Published Data on the Effect of 
Roughness on Transition from Laminar to 


Turbulent Flow 


HUGH L. DRYDEN* 


National Advisory Commitlee for Aeronautics 


ABSTRACT 


A review is presented of the published data on the effect of 
roughness, especially single roughness elements, on transition 
from laminar to turbulent flow, in which an attempt is made to 
reanalyze and correlate the available information. The re- 
analysis shows that the transition Reynolds Number of a flat 
plate with zero pressure gradient is a function of the ratio of the 
height of the roughness element to the displacement thickness of 
the boundary layer at the element, this functional relation being 
a better representation of the data than a constant critical Reyn- 
olds Number of the roughness element. Other data show that 
the effects of roughness are similar in streams of different initial 
turbulence and that a plot of the ratio of transition Reynolds 
Number of the rough plate to that for the smooth plate against 
the ratio of the height of the roughness element to displacement 
thickness of the boundary layer at the element gives good cor- 
relation of all the data for a given shape of roughness element 
when transition occurs downstream from the roughness element. 
At a certain value of the height-thickness ratio dependent on the 
stream speed, location of roughness element, and air-stream 
turbulence, the transition position reaches the element and 
remains there as the height or the stream speed is further in- 
creased. 

The paper also discusses available data on the effect of dis- 
tributed roughness on transition on a flat plate, as well as some 
of the published data on roughness effects on transition on air- 


foils. 


INTRODUCTION 


PREPARATION OF AN ARTICLE for the 


Program on 


oo THE 

Princeton Publication 
from laminar to turbulent flow, a review was made of 
the published data on the effect of roughness on transi- 
tion. The data were reanalyzed in terms of various 
concepts of the physical action of roughness in inducing 
transition and in terms of various parameters sug- 
gested by dimensional analysis in an attempt to gain 
an understanding of the fundamental aspects of the 


transition 


phenomenon. 

For the most part attention will be confined to transi- 
tion in the boundary layer of a flat plate in a stream of 
low turbulence and with zero pressure gradient. It 
is well known that in the absence of roughness transi- 
tion occurs as the result of the amplification of dis- 
turbances initially present in the fluid stream (initial 
turbulence), the effective disturbances being those 

Presented at the Aerodynamics Session, Twenty-First Annual 
Meeting, IAS, New York, January 26-29, 1953. 

* Director. 


whose frequencies lie within a specified range dependent 
on the Reynolds Number of the boundary layer. Thus, 
as indicated in the upper part of Fig. 1, transition on a 
smooth plate in a given fluid stream of mean speed V 
will occur at some distance x» from the leading edge, 
yielding a transition Reynolds Number Reo equal 
to V’xo/v, where v is the kinematic viscosity of the fluid. 
When a two-dimensional roughness element, such as a 
cylindrical wire of diameter k is placed on the plate per- 
pendicular to the fluid stream at a location x, for which 
x, is less than xo, as shown in the lower part of Fig. 1, 
transition occurs at some location x, less than x» but 
equal to or greater than x;, yielding a transition Reyn- 
olds Number for the rough plate Re, equal to Vx,/». 
The approximate course of the displacement thickness 
5* of the boundary layer is shown in Fig. 1, and the spe- 
cific value of 6* at the roughness element itself, 5*,, is 
indicated. The problem under discussion is the effect 
of single roughness elements of a given shape but of 
varying height & and varying location x, on the transi- 
tion Reynolds Number Re,. 

The published data on the effect of single roughness 
elements on transition on a flat plate are contained in 
references 1 to 6 inclusive. The data of reference 6 
were originally published in a Japanese Wartime re- 
port in March, 1945, and I am greatly indebted to 

















Prof. I. Tani for furnishing me with tables of the 
numerical data from these experiments several months 
ago. 
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Fic.2. Reynolds Number of roughness element and transition 
Reynolds Number for single roughness element on flat plate, from 
data of Tani, Hama, and Mituisi.! 


DIMENSIONAL CONSIDERATIONS 


Examination of the experimental results is facili- 
tated by the use of dimensional analysis. For the time 
being we shall consider experiments in a given fluid 
stream of constant but small turbulence on a roughness 
element of given shape. The location of transition x, 
is then a function only of the height k of the roughness 
element, its location x;, the stream speed V, and the 
kinematic viscosity v of the fluid. Since the five vari- 
ables involve only the units of length and time, dimen- 
sional reasoning tells us that the general functional re- 
lation between the five-dimensional quantities can be 
reduced to a relation between three independent non- 
dimensional parameters. As these parameters we 
might select, for example, the transition Reynolds 
Number Re, = Vx,/v, the nominal Reynolds Number 
Vk/v of the roughness element, and either the ratio 
x,/X, Or xX;/k. Many other nondimensional com- 
binations of the five quantities can be formed; any 
three independent ones may be used. Thus, rather 
than the nominal Reynolds Number of the roughness 
element based on the free-stream velocity, one might 
consider the Reynolds Number Re, based on the ve- 
locity u, in the laminar boundary layer at the height of 
the roughness element. The slope of the velocity dis- 
tribution curve at the wall at x, is given by the ex- 
pression 

(du/dy)y = 0.332 V (vx,/V)~' 
Hence, for small values of k, 
Uy, = 0.332kVy~ x, — 7 


Re, = u,k/v = 0.332 (V/v)'*k2x, 


This combination of the original variables has a more 
x,/Rk is pre- 


significant physical meaning than Vk/v. 
ferred to x,/x, as the third parameter, since the prin- 
cipal dependent variable of interest x, then appears in 
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only one of the three parameters and x;/k is determined 
solely by the geometrical characteristics of the mode] 
under investigation. Still another possible nondimen- 
sional combination of physical significance is k/],72 
(vx,/V)'? = k/6*,, the ratio of the roughness height 
to the displacement thickness of the boundary layer at 
the location of the roughness element. 

The conclusion reached by dimensional analysis js 
that the experimental results should be capable of repre- 
sentation by a relation between three independent non- 
dimensional parameters which may be selected in a 
great many ways. It may be that by a suitable choice 
of parameters the effects of one or even of two of the 
three parameters are small or absent altogether. In 
the choice of parameters we may be guided by some 
physical concept of the phenomenon or an empirical 
approach may be used. In the present analysis of the 
experimental data, the results were plotted in terms of 
several parameters suggested by physical concepts in 
the search for the simplest and most meaningful pres- 
entation. 

PHYSICAL CONCEPTS 

An early physical concept of the action of a roughness 
element in inducing transition was advanced by Schil- 
ler.” Having studied the transition in the flow behind 
a cylinder and flat plate at very low Reynolds Number 
from a laminar flow to a separated flow with vortex 
street, he suggested that a roughness element induces 
transition when the Reynolds Number of the element 
reaches a definite critical value at which vortices appear. 
The Reynolds Number used by Schiller was based on 
the friction velocity, equal to the square root of the 
ratio of the shearing stress at the surface to the density. 
It is easily shown that Schiller’s Reynolds Number is 
equal to the square root of Re,. In this theory it is 
tacitly assumed that, as a roughness element is in- 
creased in height at a given location and for a given 
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roughness element to boundary-layer displacement thickness at 
element, from data of Tani, Hama, and Mituisi.! 
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stream speed, nothing occurs until a critical height is 
reached. Then transition occurs at the roughness ele- 
ment itself. 

A second physical concept arises from discussions in 
the literature of roughness effects in pipes, in which it 
is frequently stated that the effect is dependent on the 
height of the roughness element in relation to the bound- 
ary-layer thickness.’ In addition to the intuitive 
appeal of this concept, some support is provided by the 
Tollmien-Schlichting theory of boundary-layer insta- 
bility. The roughness element imposes a disturbance 
that is added to those already present from the initial 
turbulence of the stream. Its contribution to the total 
disturbance within the critical wave-length range for 
which amplification occurs, depends on its shape and its 
height relative to the boundary-layer thickness. Its 
presence produces earlier transition, since with a 
larger disturbance less amplification is required. [If all 
disturbances are small, we might expect the roughness 
effect to be independent of other disturbances present 


regardless of their source. 


RESULTS OF TANI, HAMA, AND MITUIS1 


The early experiments of Tani, Hama, and Mituisi' 
were made on cylindrical roughness elements (wires) 
on a flat plate with zero pressure gradient in a wind 
tunnel in which Re, for the smooth plate was 2 X 10°. 
In the original paper the experiments were interpreted 
in the light of the first concept—that of a constant 
critical Reynolds number of the roughness element. 
Transition was measured by a surface Pitot tube at a 
position x, which remained the same for all of the meas- 
urements. For several values of the roughness height k 
and location x; of the roughness element, the stream 
speed was varied until transition appeared at the fixed 
location x, as judged by the sudden rise of Pitot pres- 
sure. The assumption, Re, = C, a constant, with tran- 
sition assumed to occur at the roughness element—.e., 
at x,— leads to the relation 


o~ 


Re, = uyk/v = 0.332 (V/v) ”* k?x, = 

written in reference | in the form 
0.576k/x, = C’*(Re;) 

From a logarithmic plot of Re, vs. k/x,, C * was deter- 
mined to be 13, or C = 169. A comparison of the ob- 
served data with values computed from this assumption 
of constant Re, shows differences of +40 per cent from 
the average value of Re,. 

It is known from the later experiments® that transi- 
tion does not move suddenly from its original smooth 
plate location to the roughness element as the height 
of the element is varied. Rather there is a gradual 
progression and the element may influence transition 
far downstream. Thus the experiments in reference | 
must be reanalyzed by using the actual value of x, in 


computing Re;,. 
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Fic. 4. Transition Reynolds Number and ratio of height of 


roughness element to boundary-layer displacement thickness at 
element, from data of Tani, Hama, and Mituisi' and Tani and 
Hama.‘ 


Values of Re,, Re,, and k/6*, were computed from 
Fig. 2 shows a plot of Re, vs. Re, 
The assumption 


the original data. 
and Fig. 3 is a plot of Re, vs. k/é*;. 
that Re, is constant is not a good representation of the 
data, as shown by Fig. 2. On the other hand, the as- 
sumption that Re, is a function of k/6*,. and independent 
of Re, gives a reasonably good correlation of the data. 
From these data alone, one of the three nondimensional 
parameters required by dimensional considerations 
seems to have no effect. However, the height of the 
roughness element was never increased sufficiently to 
bring transition forward to the element, since the transi- 
tion was always observed at a fixed x, greater than x,. 
Obviously, «, can never be less than x;,. 

Fortunately, Tani and Hama continued their meas- 
urements. The data in these experiments,® which were 
made available to me by Professor Tani, were analyzed 
in the same manner with the results shown in Fig. 4. 
Here it became obvious that departures from a single 
functional relation between Re, and k/6*, occurred as 
the transition position approached the position of the 
roughness element. The transition position x, does 
not move forward of the position x, of the element. 
A suitable nondimensional parameter for the roughness 
element position is x,/k. Ona plot of Re, vs. k/6*,, the 
limiting case x, = x, corresponds to the curves repre- 
senting 

Re, = 2.96(x;/k)2(Rk/6*,.)? 
since 
Re, = Vx,/v = Va,/v = 2.96 x,7/(6 


A few of these limiting curves corresponding to some 
of Tani and Hama’s values of x,/k are shown in Fig. 4. 
The experimental points begin to depart from the Re, 
vs. k/é*, relationship as x, approaches x, and then 
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roughness element to boundary-layer displacement thickness at 
element, from data of Stiiper.* 


approach the curves representing x, = x, for the perti- 
nent value of x;,/k. 

There are small systematic differences between the 
data of references 1 and 6, the cause being unknown. 
Perhaps the initial turbulence of the stream or the initial 
roughness of the plate was somewhat different in the 
two series. 

To summarize, if we select as the three nondimen- 
sional parameters the transition Reynolds Number 
Re,, the relative height of the roughness element k/6*;, 
and the roughness element location parameter x;/k, the 
data are well represented by a single curve of Re, vs. 
k/65*, for all values of k, x,, and V so long as the result- 
ing x, is greater than x;. 


RESULTS OF FAGE, LIEPMANN, STUPER, AND 
SCHERBARTH 


9 


Fage? made measurements on roughness elements 
of various shapes, including a smooth bulge, a smooth 
hollow, and a flat ridge, in an air stream with a favor- 
able pressure gradient. His data show the large in- 
fluence of the shape of the element, the flat ridge pro- 
ducing much more effect than a hollow or bulge of the 
same height. Fage’s data were, I believe, the first to 
show that a roughness element could move transition 
forward without necessarily bringing it to the element 
itself. 

Liepmann® studied conditions at a single roughness 
element in some detail, the element being a half cylinder 
with axis normal to the stream and of relatively large 
scale. Although the flow always separated from the 
surface of the roughness element, it often reattached 
to the plate with laminar-type distributions of mean 
Unfortunately, Liepmann could make meas- 


velocity. 
urements only for some 8 in. downstream, which, on 
the large scale of his experiments, corresponded only to 
a distance of the order of 50 to 150 times the roughness 
height. His 


transition Reynolds Number for the 
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smooth plate was only 500,000 due primarily to trans- 
verse contamination from the boundary layers at the 
tunnel side walls. 
than 0.7. 

Stiiper investigated cylindrical roughness elements 


His values of k/5*, were all greater 


on a plate with and without pressure gradient. His 
results are not directly comparable with those pre- 
viously discussed because he used a different technique 
for detecting transition. He recorded the variation of 
speed with time by means of a hot-wire anemometer and 
measured on the oscillogram the proportion of the total 
time of observation that the flow was turbulent. No 
matter what ratio of time of turbulent flow to time of 
laminar flow was selected as the criterion for transition, 
the value of Re, for the smooth plate decreased rapidly 
with increasing speed. The results are presented in 
small-scale diagrams difficult to read accurately. Typi- 
cal results are shown in Fig. 5, the criterion for transi- 
tion being the first appearance of turbulent bursts. 
We note by comparison with Fig. 4 that, although Re, 
for the smooth plate is in the range 0.5 & 10° to 1 X 
10° as compared to 2 X 10° for Tani’s data, a roughness 
element for which k/6*, is 0.5 reduces Re, to about 
).5 Rey as in Tani’s experiments. Apparently the 
sensitivity to roughness is not greatly reduced at lower 
values of Re) corresponding to increased initial turbu- 
lence or other source of increased disturbance. 

Later investigations with Stiiper’s equipment by 
Williams and by Willis showed the presence of sub- 
stantial leading-edge disturbances and of vibration of 
the plate, so that the interpretation of the data is some- 
what in doubt. 
with respect to the influence of pressure gradient that 
‘The sensitivity 


Nevertheless, it is of interest to note 


Stiiper drew the following conclusion: 
of the laminar boundary layer to irregularities on an 
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otherwise smooth surface was not greatly altered by a 
favorable pressure gradient in the range of Reynolds 
Number covered by the experiments and with the 
amount of turbulence present in the tunnel.” 


Measurements of the influence of flat ridges 20 mm. 
wide and varying in height from 0.1 to 1.0 mm. were 
made by Scherbarth and reported by Quick,® with the 
results shown in Fig. 6. Here, too, Rep varies with the 


speed. 
DISTRIBUTED ROUGHNESS 


Holstein® made some experiments on the effect of 
distributed roughness on a wooden plate 3.36 m. long 
and 3.6 cm. thick, with an elliptical nosepiece 18 cm. 
long. The roughness was produced by attaching 
emery paper to the surface of the plate except for the 
nosepiece, which was not roughened. The juncture 
is not described in the original paper, but in a letter 
to the author, Holstein states that “... to the best of 
my recollection the emery paper had a thickness of 
about 0.3 mm. and was the same for all degrees of rough- 
ness. The junction at the nose was well smoothed 
with plasticene, so that the step of 0.3 mm. was re- 
placed by an inclined plane of about 1 cm. length. The 
roughness measurement given in the report does not 
include the thickness of the paper.’ 

The largest value of k/6* obviously occurs at the 
beginning of the roughened portion of the plate. Hol- 
stein’s results were reanalyzed in terms of k/6*), where 
6*) is the displacement thickness of the boundary layer 
at this location. Fig. 7 shows Re, plotted against 
k/6*) omitting all points for which Re, exceeds 550,000, 
since detailed analysis shows the impossibility of sepa- 
rating effects of tunnel speed, initial turbulence, and 
initial surface waviness from the effect of added rough- 
ness for larger values of Re,. Taken at face value, the 
results indicate that the effect of distributed surface 
roughness at a given value of k/6*» in reducing Re, is 
much greater than the effect of a single roughness ele- 
ment at the same value of k/6*,. This conclusion is 
probably false because Holstein’s values of k do not 
include the paper thickness and ignore the discon- 
tinuity at the end of the smooth nose. Fig. 8 shows a 
plot of Holstein’s results on the assumption that the 
roughness height at the end of the smooth nose is the 
tabulated value of & plus 0.2 mm., the fairing being as- 
sumed to reduce effectively the discontinuity from the 
paper thickness 0.3 mm. to an effective thickness of 0.2 


nm. 


CORRELATION OF DATA 


Figs. 2, 4, 5, and 6 include data in air streams for 
which the transition Reynolds Number Reo for the 
smooth plate ranges from 0.5 X 10° to more than 2 X 
10°. Roughness effects are present throughout, hence 
for all levels of initial turbulence included. If we con- 
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Fic. 7. Transition Reynolds Number for distributed rough- 
ness and ratio of nominal roughness height to boundary-layer 
displacement thickness at leading edge of roughened portion of 
flat plate, from data of Holstein.® 
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Fic.8. Holstein’s data for distributed roughness, assuming effec- 
tive roughness height equal to nominal height plus 0.2 mm. 


sider data on cylindrical wires and flat ridges and esti- 
mate values of Rey by extrapolation in those cases 
where a direct measurement is not available, we find 
that the plot in Fig. 9 of the ratio of Re, for the rough 
plate to Re) for the smooth plate vs. k/6*), gives a rea- 
sonably good correlation of the data for which x, is 
greater than, and not too close to, x,. For a cylindrical 
or flat-ridge roughness element of height equal to one- 
fourth the displacement thickness of the boundary layer 
at the roughness element, the transition Reynolds Num- 
ber of the smooth plate is reduced by 9 per cent. This 
curve applies only when transition occurs downstream 
from the roughness element, and for a given value of 
x,/k and Rey there is a limiting value of k/6*, corre- 
sponding to the transition position reaching the rough- 
ness element in its forward travel. 


ROUGHNESS EFFECTS ON AIRFOILS 


The papers of Fage and of Tani, Hama, and Mituisi 
contain some data on roughness effects on transition 


in the boundary layers of airfoils. In general, the data 
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Fic. 9. Ratio of transition Reynolds Number of rough plate 
to that of smooth plate and ratio of height of roughness element 
to boundary-layer displacement thickness at element, single 
cylindrical (circle and square symbols) and flat-strip elements 
(triangular symbols), data from references 1, 4, 5, and 6. 


are insufficient for analysis in terms of boundary-layer 
parameters. In such flows with pressure gradient it 
is useful to consider Re, as the equivalent flat-plate 


transition Reynolds Number defined by 
Re, = 0.337(V6*/v)? 


Tani and his coworkers made a few measurements with 
sufficient data to make such a computation. At k/6*, 
= ().8, the value of Re, for the airfoil is more than twice 
that for the flat plate, which would lead one at first 
sight to believe that roughness effects are less for a flow 
with favorable pressure gradient in contradiction to 
Stiiper’s results. However, Rey for the smooth airfoil 
was not measured on the airfoil at the chordwise loca- 
tion of the rough-airfoil transition. Because of the 
favorable pressure gradient, Rey is probably much 
greater for the airfoil than for the plate, and it is with 
this value that Re, for the rough airfoil should be com- 
pared rather than the value for a smooth plate without 
pressure gradient. Fage’s paper does not include data 
on boundary-layer thickness. 

There are other experiments on roughness on airfoils 

for example, those of Loftin® in which only the critical 
Reynolds Number of the airfoil (based on airfoil chord 
and free-stream velocity) for drag rise was determined 
from wake survey measurements. Insufficient data 
are available for analysis in terms of chordwise location 
of transition and boundary-layer parameters. 

For boundary-layer flow near an airfoil, the Tollmien- 
Schlichting theory shows that the stability of the flow 
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is greatly dependent on the pressure gradient and hence 


on chordwise position. The Reynolds Number at 
which amplification of small disturbances occurs may 
vary with chordwise position by several powers of 10). 
Hence, the gradual forward travel of transition with 
height of a single roughness element, corresponding to 
variations of Re, of less than one power of 10, may not 
be present. The detailed behavior is probably mark- 
edly dependent on the airfoil pressure distribution and, 
hence, on the airfoil shape. In the absence of sufii- 
cient data we can only conclude that over certain 
ranges of airfoil Reynolds Number, the transition Reyn- 
olds Number of the boundary-layer flow is a function 
of air-stream turbulence, pressure gradient, and sur- 
face roughness and that all variables probably have im- 
portant effects regardless of the magnitudes of the 
others. More data of a more detailed nature are re- 
quired to determine whether the ratio of the equivalent 
flat-plate transition Reynolds Number of the rough air- 
foil to that for the smooth airfoil at the same chordwise 
location of transition is a unique function of k/6*, for 
a given shape of roughness element. 
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Integral Relations in the Linearized Theory 
of Wing-Body Interference 


A. H. FLAX* 


Cornell Aeronautical Laboratory, Inc. 


SUMMARY 


The reverse-flow theorem is extended to the case of a wing 
mounted on a cylindrical body of any cross-sectional shape. 
Also, two Trefftz-plane theorems relating the total lift and rolling 
moment on such a body to the corresponding values for the wings 
are derived. Based on these results, some further general re- 
lationships between integrated lifts and moments on the wing, 
the body, and the combination are obtained. It is shown that 
the lift distribution on the wing-body combination in reverse flow 
is the influence function for lift of the combination. Further, it 
is shown that the influence function for rolling moment is the lift 
distribution over the combination rolling at constant angular 
velocity in reverse flow. It is also shown that the lift induced 
on the wings by a unit body angle of attack is equal ‘o the lift 
induced on the body by a unit wing angle of attack in reverse flow 
In lifting-line theory, slender-body theory, or for any wings having 
a spanwise axis of symmetry (such as rectangular or diamond 
wings), the need for consideration of reverse flow is eliminated, 
and the results assume a particularly simple form. The results 
are valid for subsonic or supersonic flow within the limitations 


of linearized theory. 


INTRODUCTION 


_— PROBLEM OF WING-BODY INTERFERENCE in the 
linearized theory of compressible flow is one for 
which few exact solutions are available. Such approxi- 
mate methods of solution as are available involve either 
a great deal of laborious calculation or are based on as- 
sumptions that would appear to restrict the applicabil- 
ity of the results. This being the case, it appears that 
some more or less generally valid integral relationships 
in the linearized theory may be useful in analysis, even 
though they provide less than the complete solution to 
the problem. 

For incompressible flow, the wing-body combination 
of minimum induced drag has been treated by Len- 
nertz' for the case of a body that is an infinitely long 
circular cylinder. This solution has been generalized 
by Pepper? to the case of an infinitely long cylindrical 
body or multiplicity of bodies of any cross-sectional 
shape. The case of a lifting line of large span mounted 
on an infinitely long cylinder of any shape has been 
given by Multhopp,* who also indicated approximate 
corrections to allow for a noncylindrical body shape of 
finite length. Multhopp’s basic solution of the wing- 
body problem requires, however, that the body be cy- 
lindrical. For very slender wing-body combinations 
with a body of circular cross section, Spreiter® has given 
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a solution that is presumably valid both subsonically 
and supersonically. Ward‘ obtained a similar solution 
for supersonic flow which, however, agrees with Spreiter's 
only for the case of a body cylindrical in the wing region. 
Miles® has pointed out that the discrepancy between 
the Ward and Spreiter solutions in the case of an ex- 
panding body in the wing region arises from the fact 
that Ward considered quadratic terms in Bernoulli's 
equation which effect an interaction between the dis- 
placement flow around the body and the lifting flow 
around the wings. Adams and Sears’ have given a dis- 
cussion of a treatment similar to Ward’s for subsonic 
The 
case of subsonic flow about wings of low aspect ratio 


flow around slender wing-body combinations. 


mounted on an infinite circular cylinder has been an- 
alyzed by Lawrence; his analysis is valid for aspect 
ratios well above those for which slender-body theory is 
valid. For supersonic flow, a successive approxima- 
tion method for wing-body interference has been pre- 
sented and applied to a wing of high aspect ratio 
mounted on a circular body by Ferrari. Nielsen and 
Pitts’ have also given a solution of the same super- 
sonic problem in terms of a Fourier series. Except for 
the case of slender-body theory, few simple theoretically 
well-founded formulas for gross lift and moment effects 
on the wing-body combination or for the distribution 
of force and moment between wing and body have been 
presented. For the case of slender wings mounted on a 
cylindrical body, Morikawa and Puckett'! have shown 
by direct calculation that the lift distribution at con- 
stant angle of attack is the influence function for total 
lift of the wing-body combination. This result pro- 
vides a simple means of obtaining the total lift for more 
complicated cases once the lift distribution for constant 
angle of attack is known. 
sult of this type valid within the framework of the lin- 
earized subsonic or supersonic theory for wings mounted 


In the present paper, a re- 


on a cylindrical body of any shape will be derived; this 
result will prove to be valid, for instance, within the 
limits of validity of any of the theories for cylindrical 
bodies cited above. Great care must be taken, how- 
ever, in attempting to apply these results to cases in 
which the body cross section varies appreciably (viz., 
Ward, Adams-Sears). 

The results of the present paper are arrived at by 
an extension of the reverse-flow theorems!* !* for lift 
ing wings to the case of wings mounted on a cylindrical 


body. These extended theorems in combination with 


9 
-) 





484 JOURNAL OF THE 


certain Trefftz-plane theorems, relating the forces on 
the body to the forces on the wings, will provide a basis 


for the derivation of a number of integral relations in 
wing-body interference. 


THE TREFFTZ-PLANE THEOREMS 


It is convenient first to derive some general rela- 
tionships between certain integrals on the trace of the 
wing and similar ones on the body in the Trefftz plane. 
The treatment of the wing-body interference problem 
in the Trefftz plane has been discussed in some detail 
by Lennertz.! In the Trefftz plane (x = ©), the ve- 
locity potential satisfies the relationship 


o* 
oy? 


a’ 


T 3.2 oi 


0 (1) 
in the vicinity of the wing-body combination for either 
subsonic or supersonic flow. This being the case, any 
two solutions of this equation, ¢; and ¢», satisfy the 
well-known reciprocal relation 


Og» O¢1 
1s = o 1s (: 
$ oor d fp oor as 


The integrals are to be taken over the body and the 
trace of the wings in the Trefftz plane. Actually, for 
supersonic flow, the flow at large distances from the 
body does not satisfy Eq. (1), but it is always possible 
to construct an outer bounding circle in the Trefftz 
plane whose radius is small enough so that Eq. (1) is 
satisfied everywhere within it to any desired accuracy 
and yet large enough so that the contribution of the in- 
tegrals of Eq. (2) taken over the curve vanishes. It 
should be noted that, since ¢(— ©) may be taken to 
be zero, ¢( ©) is given by 


"oO 
g(o) = s oe a (3) 


According to the linearized theory, the total force 
carried by any spanwise strip of wing is given by 


| > (dp, 2 
iii nwo f y j = *) ” x) 


where the subscripts u and / refer to the upper and 
lower wing surfaces, respectively. Thus, 


f(s) i puto (du = 1) (o) 
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For the cylindrical body, the total force per unit length 
of body circumference is given similarly by 


f(s) = puog(s) (6) 


This force is, of course, directed along the normal to 
the body surface at each point. The geometrical ar. 
rangement is shown in Fig. 1. 

A relation between body lift and wing lift can now be 
obtained. Consider first the flow around the body 
alone arising from an angle of attack of the body, ‘ 
As far as the flow transverse to the body axis is con- 
cerned, this corresponds to a two-dimensional flow, 
with velocity ua at z = +, past the section of the 
body cylinder. It will be more pertinent to consider 
this flow with the uniform velocity —wa _ superim- 
posed so that the velocity at z = + © approaches 
zero. The latter flow may be characterized by the 
two-dimensional velocity potential, ¢:-(¢1. — $1;) along 
the trace of the wings is, of course, zero, since ¢ is con- 
tinuous everywhere outside the body. The boundary 
condition at the body for this flow is 


0¢)/O0n = Ua COs (N, 2) (7 


where cos (”, 2) is the cosine of the angle between the 
normal to the body and the gz direction. Along the 
wings, the value of 0¢,/0n may be denoted by g(s). 
Now consider the actual flow for the wing-body com- 
bination when the body is at zero angle of attack; the 
velocity potential for this flow is called ¢2. In this case, 
the boundary condition on the body is 0¢2/O0n = 0, while 
the lift on the wings produces a discontinuity in @» at 


the trace of the wings. Applying Eq. (2) 


i Od» ‘ Od» 
I, * On oF | a on ie 


: Og: i Od) 
2 1s on — 22) 1s (8) 
I 2 on — / (do ror) l a ¢ { 


The two integrals on the left-hand side of this equation 
are zero, since O@2/On is zero on the body and (¢;, — 


¢1,) is zero on the wings. Thus, 


— Ua J ¢2 cos (m, 3) ds = — | (do, — $2) g1(s) ds 
B w 


(9) 


Multiplying by pm and denoting g;(s)/uoa as G,(5) 
li.e., Gi(s) is the flow normal to the wing surface pro- 
duced by a unit transverse velocity of the body in the 
z direction], there results 


J f(s) cos (n, 2) ds = / f(s) Gils) ds 
B w 


The integral on the left-hand side is simply the total 
lift on the body, so that 


Lp = I f(s)Gi(s) ds 


(10) 


(11) 
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WING-BODY 


Since f(s) is simply the force per unit length in the span 
direction (not necessarily in the z direction and there- 
fore not necessarily lift) on the wings, this constitutes a 
general expression for the lift carried by a cylindrical 
body in terms of the force on the wings. 

Consideration of a simple but common case will help 
to clarify the application of Eq. (11). For wings 
mounted symmetrically on the centerline of a circular 
cylindrical body, the function G,(s) which is required is 
obtained from the well-known result for two-dimensional 
flow about a circular cylinder. For this flow, ¢; is 
given by 


¢1 = (maR?*/r) cos 6 (12) 


where 7, 9 are the usual cylindrical polar coordinates; 
§ = O along the negative z axis in Fig. 1. 


’ 1 1 Og; 
G(r) = 
UWar OO 4 90% 
‘ (13) 
R* 
Gi(s) = Gi(r) = = 
y? 
In this case, then, 
R? 
Lz = I(r) — dr (14) 
A '* 


The antisymmetrical case may also be worked out 
ina similar manner. In this case, the potential ¢; may 
be taken to correspond to the flow produced by the 
cylindrical body alone rotating with angular velocity 
w about its central axis; obviously, for a circular cyl- 
inder, this is a trivial case, since no potential flow is 


produced. For a more general class of cylinders, 


Od) 


= wr cos (s, 7) (15) 
On 


where 7 is the distance from the axis to a point on the 
surface of the cylinder and cos (s, 7) is the cosine of the 
angle between the 7 vector and a vector tangent to the 
curve s at any point. Again, ¢ is continuous across 
the trace of the wing surfaces, while 0¢;/0n across the 
trace of the wing may be represented by a function /,(s). 
The potential ¢. may again be taken to represent an 
actual flow about a rolling wing-body combination. 
In this case O¢2/On on the body will be the same as 


0¢;/0n, or 
O¢2/On = wr cos (s, 7) (16) 


Then, applying Eq. (2), 


y] ¢; wr cos (s, 7) ds = J gowr Cos (S, 7) ds + 
B B 
| (dou — G2:) m(s) ds (17) 


This may be simplified to 
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f (¢2 — gi) r cos (s, 7) ds = 
B 
a‘ hy 
| (gou — e:) — (s) ds (18) 
a w 


f 100 r cos (s,r) ds = — f(s) Hi(s) ds (19) 
JB a 


u 


or multiplying by 


where //;(s) is a function giving 0¢/0n along the wings 
for unit w and where f2(s) is equal to puo(¢d2 — di), the 
total force along any cylindrical element of the body in- 
duced by the forces on the wings. The integral on the 
left-hand side of Eq. (18) may be seen to represent the 
rolling moment of these forces, since fo(s)r cos (s, 7) is 
the element of body rolling moment due to the pres- 
sures on the elemental strip of the cylinder at s. Thus, 
the rolling moment due to pressures on the body which 
may be denoted by Jz is given in terms of the forces on 


Mz = f 4 IT\(s) ds 
Ww 


Of the Trefftz-plane theorems presented here, the 
one concerning lift is the more useful, since the rolling 
moment on the body is usually small (zero for a circular 
body) unless the body dimensions are comparable to 
the wing dimensions. It should be noted that, in the 
derivations, effects of displacement of the wing wake or 
of rolling up of the wing trailing vortex sheet have not 
been mentioned. Actually, such effects may be ap- 
proximated in many practical cases by using the esti- 
mated locations of rolled-up vortices far downstream of 
the wing in place of the (¢, — ¢,) distribution over the 
It is not, however, in line with the 


the wings by 


(20) 


trace of the wing. 
main purpose of this paper to consider such further ap- 
proximations, since these must usually be devised in the 
light of the geometry of particular configurations of 


interest. 


EXTENDED REVERSE-FLOW THEOREMS 


Consider a lifting wing mounted on an infinite cylin- 
drical body of any cross-sectional shape. For this 
particular derivation, the generators of the body cylin- 
der will be assumed to lie along the x axis and to be at 
zero angle to the air stream far ahead of the wings; 
later, the effect of body angle of attack will be intro- 
duced. The perturbation velocities in the x, y, and 2 
directions are u, v, and w, respectively, and these are 
given by 


Op Od Op 


“u= », 32S 
Ox oy Oz 
where ¢ is the velocity potential. The velocity poten- 


tial satisfies the equation 
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at 
Ox” iy" sf 


» Ob O° 


(1 — AT 0 (21) 


where J/ is the free-stream Mach Number. The pres- 
sure increment over the free-stream pressure is given to 
the linearized approximation by 


fe) 
P = —puu = —pu 4 (22) 
x 
P satisfies the equation 
or ee wr 
= I) + = 0 (23) 


Ox? Oy? 02? 
on the wing, ¢ satisfies the boundary condition 


0o/On = +w, (x, S) 


where w, (x, s) is the specified velocity normal to the 
wing surface which is proportional to the local angle of 
attack, s being an arbitrary axis in the plane of the 
wing which is perpendicular to the x axis. The two signs 
correspond to the upper and lower surface values. On 
the body, 0¢/0n = 0, m being the direction normal to 
the body surface. In subsonic flow, ¢ is a continuous 
function everywhere except on the wing and in the wake 
trailing behind the wing tox = ©. Across the wing 
and wake, ¢ in general undergoes a jump in value. In 
supersonic flow, @ may be discontinuous not only across 
the wake but also across characteristic surfaces (Mach 
waves), but as in the case of reverse-flow theorems for 
wings alone,!” !* it is found that these latter surfaces of 
discontinuity make no contribution to the final results. 
Demonstration of this fact for the wing-body case 
differs in no essential way from the demonstration for 
wing alone and will therefore not be considered further 
in this paper. Also, as in the analysis for wings 
alone,!* !* consideration must be given to the singular 
behavior of P near a subsonic leading edge when it goes 
to infinity as 1/+/d, d being the distance from the edge. 
P must be continuous across the wake, since the wake 


can support no pressure difference, but it is in general 


SCIENCES SULT, 19862 
discontinuous across the wing surface and body. Jp 
order that P go smoothly into its value in the wake, 
there must also be no discontinuity in pressure across 
the wing surface as the wing trailing edge is approached, 
The flow is assumed irrotational so that 


Ou Ow Ou Ov 


= and mae (24) 
Oz Ox Oy Ox 


For a reversal of flow direction, the potential and 
pressure, @ and P, respectively, satisfy the same differ- 
ential equations as ¢ and P, but the conditions at lead- 
ing and trailing edges are interchanged. It is con- 
venient to keep the directions of positive x, y, 2 and 
positive perturbation velocity %, 3, @ the same in the 
reverse-flow case as in the direct-flow case; thus the 
perturbation pressure P is given by 


p = 


— puou (Zo 
The boundary conditions are 


Og 


on the wing and 
0¢/On = 0 


on the body. 

To’ make the discussion somewhat more specific, the 
configuration shown in Fig. 2 will now be considered. 
This consists of a body with two wings, not necessarily 
in the same plane. If the wings are not in the same 
plane, however, it must be noted that two different s 
axes must be chosen, one in the plane of each wing; 
this does not affect the procedure to be given.  Fur- 
ther, it will be clear that the reverse-flow theorem to be 
derived will not be altered if any number of additional 
wings be mounted on the body; it will, of course, be 
necessary in the case of such multiple-wing arrange- 
ments to consider the surface integrals over wings and 
wing wakes to be taken over all wings and wakes. 
Further discussion in this paper will, however, be limited 
to the configuration of Fig. 2. 

Consider now the volume integral 


. . ms = ( 4 . =) 
Ih = {P(V°o¢ — M _ 
FITC. 
; yep ve =) lv (26 
o( Ox? ]f - — 


where v, the volume of integration, is taken to consist 
of all space outside the wing-body combination. Eqs. 
(21) and (23) are satisfied by @ and 6 and P and P, 
respectively, everywhere in this space except in the 
respective wakes. Further, as x approaches + © in 
the wake of the direct flow, ¢ and P go to zero; while, 
as x approaches — © in the wake of the reverse flow, 
@ and P go to zero. At subsonic leading edges where 
the derivatives of ¢ are singular, P is zero because of 
the Kutta condition in the reverse flow (since the trail- 








ing ed 
similat 
tial in’ 
in the 
wings 
and P 
integr: 
Green 


lJ 


Here 1 
of the 
sector 


I, 


and b 


where 
norm: 
Thus, 


or 


Ip 


II 


The | 


(n, xX) 


By E 


so th 


Fron 


givin 


is 





- In 
wake, 
LCTOSs 
Ched, 


(24) 


| and 
iffer- 
lead- 

con- 
and 
1 the 
> the 


the 
red. 
rily 
une 
it s 
ng; 
‘ur- 
» be 
nal 
be 
ge- 
nd 
es. 


‘ed 











WING-BODY 


ing edge in reverse flow is then also subsonic) and, 
similarly, foréand P. Thus, Green’s theorem and par- 
tial integrations may be applied to the volume integral 
in the two spaces bounded by the planes containing the 
wings and by the body. Since ¢ and ¢ satisfy Eq. (21) 
and P and P satisfy Eq. (23) in these spaces, the volume 
integral is equal to zero over each of these spaces. By 


Green’s theorem, 


[If (PV°6 — oV?2P) dv 
Sf? P Pda — Sf fet 27) 
on * 


Here the surface of integration A is the inner boundary 
of the fluid consisting of the planes of the wings and the 
sector of the body cut out by these planes. Also, 


2 2) 
I, = —M? P —_ lv 
: SSS tee Ox? ® ax? 
af) 
= —\/[? dv (28) 
ri ie E Ox Ge Ox ia \ 


and by partial integration 


eo ae oP 
I, = M? / / (P ° —@¢ ) cos (n, x) dA (29) 
: ae a Ox Ox 


where cos (7, «) is the cosine of the angle between the 


normal to the bounding surface and the x direction. 


Thus, 


L=if,+.=90 


a - 7 op 
= / P Y us aa / / @ 
“ SA On ‘ad JA oO” 
‘Sc t.8 oP 
M? Jf (P , —@ ) cos (n,x) dA =0 (30) 
r A Ox Ox 


The last integral in this equation vanishes because cos 
(n, x) is zero over the entire bounding surface. Thus, 


> 7 m Pa . 7 oP 
/ / r . 1A / / o dA =0 (31) 
J Ja Om J Ja On 


By Eq. (2 
pity (Ou /On) 


or 


dA — 


OP /on = 


so that 


fp as [ it ae 
Z = Uy Fa = 
la) an ii a” 


From the irrotationality conditions, Eq. (24), 
Ou/On = (0/Ox) (0¢6/On) 


giving 


> 7 _ de —" 2 (28) 
¥ 1A — 0 dA =0 (33 
J J On / | ened Ox \On oo) 
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Integrating the second integral by parts with respect 
to x gives 


x 


~ Oo Og 
P 1A — 0 1s 
JS at, ¢ puy | el. | ds + 
Od OG 
lf ply 4 : dA 0 
J Ox On 


The second integral vanishes since ¢ = 0 at x = 
—o and 0¢6/0n = 0 at x = +o. Further substi- 
tuting from Eq. (22), 


<2 ° f 
ff. P * aA _ Jf r’ P aA = 0 
7 A, On A, On 


The area of integration is now only over the plane of 
the wings outside the body, since 0¢/0n and 0¢/0n are 
Eq. (35) is applicable to each of 


(34) 


(35) 


zero over the body. 
the two sub spaces into which the plane of the wings 
and the body surface divide the entire space. The sur- 
faces of integration involved are the planes of the wings 
in each case. The velocity normal to the plane of the 
wings is continuous so that 


(O0o/On), = —(0d/On) 


where the plus and minus subscripts denote upper and 
lower sides of the plane of the wings. Adding Eq. (35) 
for the two sides of the plane 


aa = ra) 
[/ (P, — P ( *) dA — 
x wa On/ , 
/ / ) (5) dA =0 
a a On/ . 


The area of integration has now been reduced to the 
wing surfaces proper, since the pressure difference across 
the wing planes (P, — P_) is zero except over the 
wings. The pressure difference (P, — P_) 
denoted by z, and w, = 0¢/O0n may be chosen by con- 
vention as being the value on one or the other side of 


may be 


the wing planes, so that 


/ / 7 wv, dA = JJ rw, dA 


This is the extended reversal flow theorem for wings 


(36) 


mounted on cylindrical bodies. 


Bopy AT ANGLE OF ATTACK 


When the body is at an angle of attack, ag, instead of 
being aligned with the flow as postulated in the pre- 
the flow around the wing-body com- 
two coim- 


vious sections, 
may be considered to consist of 
the first a uniform cross flow transverse to the 


bination 
ponents: 
body corresponding to a velocity maz at = 
which is superimposed a uniform flow of velocity — a 
at s = —o; the second a flow around the wing-body 
combination, with the body at zero angle of attack with 
For the second flow, however, 


—o,on 


velocity wm atx = —o. 
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the boundary condition over the wings must be taken 


to be 
0¢/On = w,(x, Ss) + uaz [Gi(S)] (37) 


where [G,(s)] is the dimensionless velocity along the 
wings produced by the first flow (i.e., the cross flow). 
Applying the reverse-flow theorem to a wing-body com- 
bination with the body at an angle of attack, there re- 
sults 


Sf T W,(X, S) dA a a j J T uoap|G1(s)] dA = 
[ mw @,(x, s) dA + Sf T Uya&p(Gi(s)] dA (38) 


Denoting the local geometrical angle of attack of the 


wing by a,,(x, s), 
(39) 


Wy (%; 5S) = uM gl X; S) 


so that 


JS r Ay(x, s) dA + JS g ap(Gi(s)]| dA = 
Jf wT Qy (x, s) dA +f wt &p(Gi(s)] dA (40) 


From Eq. (11) it may be noted that 


Ls = / / r [G,(s)] iA 
ai nr |G,(s)] ia\ 


where Lz and Lz are the total body lifts in direct and 
Thus, 


(41) 


™~s 
\ 
= 

II 


reverse flow, respectively. 


/ i T Ay\X, S) dA + Lrap = 
/ / © QutX, S$) dA + Lear ( {2) 


This equation may be immediately applied to obtain 
two important results in the linearized theory of wing- 
body interference. 


(A) Influence Function for Lift 


Let &, = cos 6, and & = | and let a,(x, s) be an 
arbitrary function of x and s and ag be an arbitrary body 
angle of attack. @, is the angle that the plane of a wing 
makes with the plane perpendicular to the plane of ro- 
tation of the body (the x-y plane), as shown in Fig. 1. 
It should be noted that, if the wings are held at zero 
incidence to the generators of the body while the body 
is pitched to unit angle of attack, cos 0, will be the wing 
Then, applying Eq. (42), 


angle of attack. 
|e x us Clk, S) dA + Leap = 


cos a ff wmdA+Lp, (43) 
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Since 6, is also the angle between the normal to the 
plane of the wing and the plane of the lift vector, the 
first term on the right-hand side is simply the contriby. 
tion of the wing to the lift. Of course, if the two wings 
lie in different planes and have different 0,,’s, this term 
must be considered to be the sum of two such terms. 
Finally, then, 


| o- + Le = Pes w Q(X, S) dA + Leap (44) 


which gives the total lift on the wing-body combination 
with any distribution of wing angle of attack and any 
body angle in terms of the lift distribution of the same 
wing-body combination at unit angle of attack (with 
zero wing incidence) in reverse flow as an influence 
function. It should be noted that, if @, is not zero, 
L» refers to the component of force on the wing in the 
lift direction, not to the total force on the wing. 
wing symmetrical about a spanwise axis, such as a rec- 


For a 


tangular wing, or within the scope of lifting-line theory 
or slender-body theory, the introduction of the reverse 
flow is unnecessary. This result may be regarded as an 
extension to wing-body combinations of the generalized 
Munk formula for lift previously given for wings alone." 
The result of Morikawa and Puckett!' is a special case 
of Eq. (44) restricted to very slender wings mounted on 
circular bodies. 


(B) Reciprocity Between Induction of Lift on Wing and 
Body 


Let &, = Oand &@g = 1, while a, = cos 0, and ag = 
0. The reverse flow then corresponds to the case in 
which the body is at unit angle of attack, while the 
wings are adjusted to a negative constant incidence 
so that their geometric angles of attack are zero. The 
direct flow corresponds to the body being at zero angle 
of attack, while the wings are at angle of incidence cos 
6,. (It must be understood here that all angles of in- 
cidence and angles of attack must be small enough to 
conform to the limitations of the linearized theory. 
The unit angle of attack referred to is simply for mathe- 
matical convenience as are the angles cos @,.) Apply- 
ing Eq. (44), 

= Lp (45 
Thus, the lift induced on the wings by a unit body angle 
of attack with the wings at zero angle of attack is equal 
to the lift induced on the body (at zero angle of attack 
by the wings at angle of attack cos @,. It should be 
noted here that a wing angle cos 6,, corresponds to the 
wing angle of attack which would be produced geo- 
metrically (with zero wing incidence) by a unit body 
angle of attack. 

Correlations of these results with extensive experi- 
ments on supersonic wing-body interference have been 
given by Cramer.'' In view of the errors introduced by 
finite body length and the gaps at the wing-body junc- 
ture for the wing incidence cases, the agreement ob- 
tained was fairly good. 
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WING-BODY 


ANTISYMMETRICAL LOADINGS 


For antisymmetrical loadings, the total rolling mo- 
ment replaces the total lift as the item of interest. For 
the special case of a circular body, the pressures on the 
body produce no rolling moment, so that only the con- 
tribution of the wings is of interest. In this case, Eq. 
(36) may be applied directly to obtain the desired re- 
sult. For noncircular bodies, however, the body may 
contribute a rolling moment, and the Trefftz-plane 
theorem for the antisymmetrical case must be applied 
in combination with the reverse-flow theorem to ob- 
tain the influence function for rolling moment. Con- 
sider the case of a wing-body combination rolling with 
angular velocity w about an axis in the body. As in the 
symmetrical case, this flow may be considered to be 
made up of two constituent flows. The first of these is 
the flow with potential ¢, corresponding to the flow 
about the infinite rotating cylinder without wings. 
The cross-flow velocities corresponding to this flow 
satisfy the boundary condition everywhere on the sur- 
face of the cylinder. On this must be superimposed 
a second flow satisfying the boundary conditions on the 
wings while giving rise to no additional flow transverse 
to the surface of the body. The wing boundary condi- 
tion in this latter case must be, for unit angular velocity, 


O¢2/On = r cos (s, 7) + H(s) (46) 


where //,(s) is the velocity normal to the wing plane 
produced by ¢, for unit w and cos (s, 7) is the cosine of 
the angle between the radius vector from the axis of 
rotation to a point on the wing and the plane of the 


wing. Also, on the body, 


O¢2/O0n = 0 
Thus, the reverse-flow theorem, Eq. (36), may be ap- 
plied to the flow represented by ¢2. Now, let @,(x, 5) = 
cos (s, r) + H,(s) and let w,(x, s) be any other dis- 
tribution of specified normal velocity over the wings. 
Then, by Eq. (36), 


JIS wT W,(X, S) iA = 
/ / a [rcos (s,r) + A, (s)] dA (47) 


According to Eq. (20), the rolling moment on the body, 


Mp, is given by 


w= ff nw I1,(s) dA 


while it is evident that the rolling moment on the wings, 
M,, is given by 


M, = / / mr COS (S, 7) dA 


Thus, the total rolling moment on the wing-body com- 
bination is given by 


(48) 


(49) 
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M, + Mz = | | T W,(x, Ss) dA (50) 
wheie z is the pressure distribution over the wings in re- 
verse flow due to a unit rolling velocity of the wing-body 
Clearly, z is the influence function for 
Eq. (50) can again be con- 


combination. 
total rolling moment. 
sidered to be a generalization of a similar result for the 


wing alone.!® 
CONCLUSION 


A number of integral relations for wings mounted on 
cylindrical bodies of any shape have been obtained on 
the basis of the extended reverse-flow theorem in com- 
bination with some Trefftz-plane theorems. Obvi- 
ously, the extended reverse-flow theorem has further 
applications similar to those discussed for wings alone, '* 
but these will be apparent from the corresponding work 
for wings. Further, the adjoint variational principle 
given for wings alone!® may readily be extended to wing- 
body combinations through the use of the extended re- 
verse-flow theorem; the demonstration of this fact 
follows in such close parallelism to the wing-alone case 
that it is not deemed necessary to present the detailed 
argument in the present paper. The Trefftz-plane 
theorems relating to the lift and rolling moment due 
to the body to those due to the wings have been found 
to be of considerable value in setting up and estab- 
lishing the limits of error of various approximate meth- 
ods for analysis of wing-body interference, and in the 
interpretation of experimental results. The details of 
such methods are strongly dependent on the geometry 
of the configuration being studied and, because of this, 
are also not included in the present paper. 
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Asymptotic Forms of Shock Waves in 
Flows over Symmetrical Bodies at Mach | 


DAVID T. BARISH* ann GOTTFRIED GUDERLEY* 
Wright dir Development Center 


SUMMARY 


The forms of shock waves in flows at Mach 1 at large distances 
from a symmetrical body are considered for planar and for axi- 
symmetric flows. It has previously been shown that for the 
flow upstream of such shocks similarity solutions exist. It is 
found that the assumption of similarity solutions downstream of 
the shocks is consistent with the shock relations in the transonic 
approximation. These facts are used to compute the flow fields 
at great distances from the bodies. 


(1) INTRODUCTION 


——— OF MIXED FLOWS over the upper half of a 
symmetrical body with various free-stream veloci 
ties close to sonic are shown in Fig. 1. In each example 
a supersonic region extends from the body out into the 
flow field. At low subsonic Mach Numbers with 
mixed flows, the supersonic region may terminate on 
the body; at higher subsonic Mach Numbers it ter- 
minates somewhere behind the body; 
sonic Mach Numbers it will not terminate at all but 
In general, it 


and at super- 


will extend downstream to infinity. 
can be expected that the return from supersonic ve- 
locities to subsonic velocities will be brought about 
If the supersonic flow terminates on the 
At low supersonic 


by a shock. 
body, exceptions may be possible. 
velocities, there is a detached shock in front of the 
body. 

At Mach Number 1 (Fig. 
contain neither the shock in front of the body, which is 


lc), the flow pattern will 


characteristic for supersonic Mach Numbers, nor the 
shock downstream of the body, which at high subsonic 
Mach Numbers region. 
However, another type of shock must be expected in 
flow patterns at a Mach Number 1 and close to 1 for 


terminates the supersonic 


the following reason: 

For convenience, consider the upper half of a planar 
symmetrical flow. The structure of the supersonic 
field is such that the Mach waves traveling toward the 
sonic line are rarefaction waves. The last of these 
Waves to reach the sonic line has been called the ‘‘limit 
ing’ Mach wave. Those waves that reach the sonic 
line will be reflected at the sonic line as compression 
waves. 

Now assume that the body shape is such that down- 
stream of the starting point of the limiting Mach wave 


no waves are reflected from the body. Then down 


Received August 9, 1952. 
* Flight Research Laboratory. 


stream of the limiting Mach wave the flow field contains 
only one family of Mach waves—i.e., those Mach waves 
that start at the sonic line. In this field the velocity 
vector is the same all along one characteristic, and its 
value can be found at the point of the intersection of 
the limiting Mach wave with the Mach wave con- 
sidered. Since the velocity vector at the limiting Mach 
wave tends to become sonic and horizontal as one goes 
out further and further, such a flow will tend to become 
parallel and sonic. This body, however, will have an 
ever increasing thickness in the downstream direction. 
For an actual body of finite length, it can be seen that 
waves will be emitted; however, the sum of the rarefac- 
tion and of the compression waves downstream of the 
limiting Mach wave is zero. (One can regard the axis 
of symmetry as a rigid wall at which waves are re- 
flected.) The waves reflected from the body will de- 
pend upon its shape. First, one will have rarefaction 
waves that turn the stream lines back toward the axis 
of symmetry. Further downstream, one will either 
have compression waves that may start at a portion of 
the body with reduced curvature or one will have a 
shock that starts at the trailing edge. Since the waves 
reflected at the line of symmetry are all compression 
waves, the sum of the waves that start at the body 
will correspond to a rarefaction. The compression 
waves that start at and behind the body will coalesce 
toformashock. In the immediate vicinity of the body, 
the shape of the shock depends upon the body shape, 
but at a distance from the body the shock is essentially 
determined by the compression waves that emanate 
at the sonic line at a large distance from the body and 
then are reflected at the axis of symmetry. The influ- 
ence of the body shape on the distribution of these 
waves is small. Accordingly, the form of the shock at 
a great distance will not show this influence either. In 
the present paper, the shape of this portion of the shock 


will be determined. 


In a flow pattern with a free-stream Mach Number 
1, the sonic line extends from the body to infinity. 
Among the Mach waves that start at the body, there is 
the limiting Mach wave mentioned above. Perturba- 
tions introduced downstream of the limiting Mach wave 
will not influence the flow field upstream of it. The 
flow field upstream of the limiting Mach wave has been 
determined previously. It has been shown that, at 
large distances from the body, the velocity potential ¢ 
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FORMS OF 


which describes the deviations from a sonic parallel 
flow has similarity solutions of the form 


—e 


3n-2 ¢ 
¢=y" “f(xy ") 


where x and y are the physical coordinates. 

Now consider the possibility of describing the flow 
behind the limiting Mach wave at a great distance also 
The two flow patterns would 


by a similarity solution. 
Further- 


have to match along the limiting Mach wave. 
more, at the line of symmetry downstream of the body 
the velocity component normal to the line of symmetry 
would be zero. A detailed study of such solutions has 
shown that a shock-free flow of this kind does not exist, 
which is not astonishing in the light of the previous 
description of the flow pattern. 

The basic idea of the present paper is the following: 
If a shock is assumed to occur along a line x-y ” 
constant and if one carries out the same simplifications 
in the shock conditions as in the flow differential equa- 
tion, then the flow after the shock is again representable 
by a similarity solution. When such a shock is ad- 
mitted, a flow can actually be found which fulfills all 


physical conditions of the problem. 


(2) SIMILARITY SOLUTIONS 


In the general approach used here, the analysis for the 
planar and axisymmetric cases are combined. «x is the co- 
ordinate in the free-stream direction; y is the coordinate 
normal to it; y is the ratio of the specific heats; a* is the 
sonic velocity, and A is a constant that is 0 for planar 
flows, and | for axisymmetric flows. The velocity po- 


tential @ may be written in the form 


a*(x + ¢) 


() = 
i.e., ¢ describes the deviation of the flow from a sonic 
parallel flow. The equation for ¢ can be simplified by 
the considerations of transonic similarity, and one ob- 


tains 
(y + 1) ¢rGrr — Sw — Al(y,/y) = 0 (1) 
The solutions considered will have the form 
g = y"-? £(¢) (2a) 
where 
f=(74+1)"%xy ° (2b) 


For f one obtains the following ordinary differential 
equation: 


fine? — f') — ngf'(dn — 5 + A) + 


f(38n — 3 + A) (8n — 2) =0 (3) 


The primes denote derivatives with respect to ¢. 

It has been determined that for the planar case n = 
4+ 5and in the axisymmetric case n = 4/7. The above 
differential equation must be integrated numerically. 
It will be noted that, for the limiting Mach wave, the 
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coefficient of the term f” vanishes—i.e., f’ = n°g*. ¢, 
for the limiting Mach wave, may be chosen to be 
equal to 1; then, f(¢) can be expressed by the series 


fg) = > Axe — 1)" (4) 
A=0 


where the coefficients are: 


In the Planar Case In the Axisymmetric Case 


Ay = 2.1333 Ap = —2.61224 
A, = 0.6400 A, = 0.32653 
Ay = 1.1200 A, = 0.73469 
A; = 0.25315 A; = 0.20408 
A, = 0.02115 1, = 0.00512 
A; = 0.39088 1, = —0.001408 


For large positive values of ¢, the function, f, can be 


expressed in the form 


ri ‘se ‘5. \ a %e\(a—3)Ar - 
f(Ce) = (CE)* & BCS) (5) 
A=0 


where C is an arbitrary constant. Substituting this 


expression into the differential Eq. (3), one obtains for 


for an axisymmetric flow, a = 


a planar flow a = ! 9; 
l 


A change of the value of Bo is equivalent to a change 
of C. If B is chosen equal to 1, the coefficients become 


For Planar Flows For Axisymmetric Flows 


By = +1 By = +1 

B, = —0.06250 B, = —0.09375 
B, = +=0.018229 B, = 0.07617 
B; = —0.009766 B; = 0.09219 
By, = 0.006744 By, = 0.03490 
B, = —0.005372 


With the upper signs, one obtains supersonic veloci- 
ties; with the lower signs, subsonic velocities. These 
series developments will be used in order to obtain 


starting points for the numerical integrations of Eq. (3). 


(3) THE s, t PLANE 


Eq. (3) possesses a group property that can be ex- 
pressed in the following manner: If fi) = g(f) is a 
solution of the differential equation, then fi2, = C*g(Cf) 
Here, C is an arbitrary constant 
corresponding to C in Eq. (5). Solutions that differ 
only by the value of C are related to each other by a 
By the 


is also a solution. 


scale transformation in the physical plane. 
transformation 

(6) 
the constant C drops out; therefore, all solutions in 
the physical plane which differ only by the scale will 
coincide if one expressed them by s and ¢. The intro 
duction of s and ¢ reduces the order of the differential 


equation, and one obtains 


dt 

ds 

12 + (3n2? — 5n + An) — s(3n — 2) (38n — 3 + A) 
(n* — t) (t — 3s) 
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FORMS OF SHOCK WAVE 


The representations in the s, ¢ plane (Fig. 5) give a use- 
jal survey of the properties of the solutions. The fol- 
lowing facts are of importance: 

In the finite part of the plane, there are three singu- 
lar points—A, B, and C. They have the following 


coordinates: 


A B ay 
‘ 0 n3(5n — 5 + A) 
= 2\/2 ‘ (3 — A)/9 
(3n — 2) (8n —3 + A) 
t 0 n? (3 — A)/3 


At the singular point A, the value of ¢ tends to in- 
This point, therefore, corresponds to the hori- 
A is a nodal point. 


finity. 
zontal axis in the physical plane. 
The integral curve in the s, ¢ plane which is obtained 
irom Eq. (5) is the only one that gives, along the x-axis, 
no vertical velocity components. Depending upon 
whether one takes the portion of the curve in the s, / 
plane which lies above or below the x axis, one has 
either supersonic or subsonic velocities. 

The singular point B lies at a finite value of ¢; it cor- 
responds to the limiting Mach wave. Point B is a 
saddle point; hence, only two integral curves enter 
this point. One of them goes to point A; the other one 
goes to point C. One might expect that the curve BA 
in Fig. 5 represents the solution that prescribes the flow 
field downstream of the last Mach wave. However, it 
is found that the integral curve BA enters point A in 
such a manner that the vertical velocity components 
along the axis in the physical plane are not zero. The 
curve BA represents the aforementioned shock-free 
flow behind the limiting Mach wave; it is the flow over 
a suitably shaped half-body. 

As the point C is approached along an integral curve, 
the values of ¢ go to infinity—1.e., this point also corre- 
Since s and / are finite, this means 
Accordingly, the 


sponds to the x axis. 
that f and f’ for point C are infinite. 
velocities along the x axis in the physical plane are in- 
finite. How can such infinite velocities arise? The 
supersonic flow downstream of the limiting Mach wave 
can be constructed by means of the method of char- 
The Mach waves that cross the limiting 
Therefore, these in- 


acteristics. 
Mach wave are predetermined. 
finite velocities can be obtained only by rarefaction 
waves that start at the body—i.e., in our case, at the 
origin of the x, In other words, the curve 
BC corresponds to a Meyer expansion, which occurs 
downstream of the limiting Mach wave. That 
can reach velocities that are infinite, and that the Mach 
wave for which these infinite velocities are obtained is 
horizontal, comes from the fact that the simplified flow 


y system. 


one 


equation, Eq. (1), has been used. 

The y-axis in the physical plane maps to infinity in 
the s, ¢ plane. The integral curve in the s, ¢ plane 
which describes the flow upstream of the body starts at 
point A. i.e., to the point 
that corresponds to the y axis. 
out the actual integration in the x, y plane; then the 


It then goes to infinity 
(It is simpler to carry 
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transition through the y axis does not offer any com- 
plication.) With the proper choice of the value of n, 
the curve will then go to point B. The analytical 
continuation of this integral curve beyond B 


is the integral 


1.2., 
downstream of the limiting Mach wave 
curve BC. 

Although it is plausible that the flow downstream of 
the limiting Mach wave is the analytical continuation 
of the flow upstream, this is not a necessary physical 
requirement, since it is possible that along the limiting 
Mach wave a singularity propagates. Therefore, it 
would be justified if one discusses the integral curve BA 
as a possible solution for the flow after the limiting 
Mach wave. However, it was mentioned previously 
that this curve does not give zero velocities along the 
x axis. Also, the curve BC as a whole is unsatisfactory, 
since at point C it gives infinite velocities along the 
v axis. 

The following analysis discusses whether it is possible 
for the actual solution to follow first the integral curve 
BC and then to jump, by means of a shock, to another 
integral curve that gives zero velocities along the x axis. 


(4) SHock RELATIONS 


The following additional symbols will be used: wis the 
modulus of the velocity vector and e is the inclination of 
The subscripts J and // de- 
note conditions before after the shock, 
The subscript NV denotes the components of 
It 
sumed that the shock occurs along a line ¢ = ¢y. 

The velocity vector in front of the shock is described 


the normal to the shock. 
and respec- 
tively. 
the velocity vector normal to the shock. 


is as- 


by 
G2, = By"*~* (Sa) 
with B = f,(¢)) and 
i = Dy?" 3 
with 
D = (8n — 2)f1r(o) — no&fr'(&o) (Sb) 


From the definition of — [Eq. (2b)], it follows that the 
slope of a normal to the shock is given by 


tan e = Ey"! (9) 


with £ The slope of the velocity vector is 
given by ¢,,. 
that, if one travels along the shock (i.e., along ¢ = fo) 
in the direction in which one approaches the sonic 
state (i.e., for values of nm > 1 toward small values of 


— nfo. 
A comparison of Eqs. (9) and (Sb) shows 


y for values of nm < 1 toward large values of y), the 
angle between the velocity vector and the normal to 
the shock is given essentially by e. 

In front of the shock, the component normal to the 
shock is therefore given by 


Wy Wy COS € 
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FIG. 4 FUNCTIONS FOR AXISYMMETRIC FLOW 
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In the lowest order approximation, it is found that one 


obtains, from Eq. (9), 


9 /y) 4,22 — 1 
cos e = 1 — (E?/2)y"" 


Wry - on —2 B | i 
a” ere é z 


Bernoulli’s equation in the form 


Thus, 


ay” Y + l oe l Wy 

a*? 2 3 wt 
vields 

ay" Qn > 

iis ae” 

gg” 


The normal component behind the shock is expressed in 


terms of the condition in front of it by 


WIN = 2 ay" ay" Y 


a y+ la*® wy, y+1a,* 


| Wry 


(10) 


Introducing the quantities a," and wy,, from the previ- 
ous equations, one obtains 


WIT, vy 


=]-— 


- : B(y + 1) + - (y — 3) 


Therefore, the velocity change normal to the shock 
and, since the tangential velocity components do not 
change as one passes through a shock, also the entire 


velocity change in the shock is given by 


Wiy — Wiry a Dye" » (3 a fk? ) 
a* ie y¥+1 


‘ro is one can determine the change of the x. an 
From this on n determine the change of th nd 


ycomponents of the velocity 


(Wry — Wry) ‘ “ 
J —_ you - 
nun, = O11, — * CoS <€¢ = 5 x 
a 


(Wr, 


Si ss E+ 
vy “| D—2V2E ( _ ) 
: y+ 1 


If the flow after the shock has also similarity solutions 


of the form 
* fil) 
then, 
en: = ¥"" fr’ 
* (3a — 2) fir — nv fir’ 
Equating in the last formulas the values of ¢g,, and ¢,,, 


one obtains, by the transformation to the variables s 


and /, the following shock conditions: 


OF SHOCK WAVES IN FLOWS AT 


MACH 1 $97 


thy = 2n- — ty (lla) 


(1l1b) 


In finding the shock locations, the use of the s, ¢ plane is 
particularly advantageous, since, in the present prob- 
lem, the scale of the postshock flow is unknown. In 
the s, ¢ plane, this scale factor does not appear. 

For other examples of two-dimensional flow, one 
may be interested in carrying out the entire analysis in 
the hodograph plane. 


(5) Tue Actua, FLow PATTERN 


The flow pattern with the shock is found by the fol 
lowing procedure: 

(a) One determines the preshock flow downstream of 
the limiting Mach wave. This is done by a numerical 
integration using Eq. (4) in order to obtain a starting 
The result is then transformed into the s, ¢ 
plane. This is the path from B toward C in Fig. 5. 

(b) For these values of s and ¢, the curve for the con- 
ditions after the shock is found by means of Eq. (11). 
This curve is the locus of all possible conditions behind 


point. 


shocks that occur along surfaces of constant ¢. 

(c) The postshock integral curve is obtained, using 
Eq. (5) to determine a starting point and then pro 
ceeding numerically. The result is represented in the 
s, t plane by the line emanating from A toward B. 

(d) The intersection of the curves found under 2 and 
3 establishes the shock location and determines, in the 
flow upstream of the shock, the value of ¢) and, in the 
flow behind the shock, the value of (cf). 

(e) From these values the quantity C is determined. 
The following results were obtained: 

s t Cc 
0.5185 0.330 1.579 2.018 


Planar flow 
0.0610 0.032 1009 2 240) 


Axisymmetric flow 


(6) THE SHAPE OF THE STREAM LINES 


Since, in the simplified potential Eq. (1), only small 
deviations from a rectilinear flow are considered, the 
deviational ordinates, ¥, of the stream lines from the 
stream lines of the basic sonic flow can be found by 
integrating the stream-line slope along a line vy = con- 
stant = VW. 

This integral can be expressed in terms of the f and f’ 
by the following considerations: 

7lanar case: Ina region S with the contour C, one 
has, because of the flow differential Eq. (1), 


Jf [(y + 1) Grrr — Gy, | dx dy 0 


By Green's theorem, one obtains 


J (y + 1) dy + [eva = (0) 
J o at 


The direction of travel around this contour is counter- 


clockwise. If one chooses in particular a region bound 
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by the lines y = 0, vy = Yo, x = — ©, and ¢) = constant, 
all terms drop out except 


*x0 7K o,* 1x 
yg, dx = Ec + 1) + gy | dy 
a 0 Zz dy 


hence, for n = 4/5, 


ie wee Se. - 
y= / g, dx = Sy” E = = (f — 2éf | 
= yA a0 


Axisymmetric case: By an analogous procedure, 


which starts with the integral 


J / (y + l)ergrr — Sw — ¥y/y] dx dy = 0 
s 


one obtains the expression 


~ Se 7 - rr S & 4 .ef/ 
+, = y dx = > 9 = 19 Ef = 2éf ) 


The forms of the stream lines are shown in Fig. 2; actu- 
ally, the transonic approximations are valid only for 
larger values of x and y (smaller perturbations) than 
those shown in Fig. 2. The values of the deviational 
ordinates and the pressures are given in Table 1 and 
are shown in Figs. 3 and 4. 


(7) DISCUSSION OF THE RESULTS 


Generally speaking, the computations show that, at 


a sufficient distance from a body, not only the flow in 


INTEGRAL CURVES IN THE S,T PLANE 


front of the limiting Mach wave but also downstream 
is described by similarity solutions. This implies that 
the flow field has a form that is independent of the body 
shape and corresponds to the discussions about the 
form of the shock made in the introduction. 

It will be noted that in the planar case the deviational 
ordinate of the stream line is proportional to y ‘—1.e., 
it tends to infinity as y tends to infinity. However, 
the stream-line slopes go to zero as one goes further out. 
This fact is, to some extent, in contrast with the usual 
concept of the properties of the flow pattern. It shows 
that the displacement of the stream lines at a large dis- 
tance from the body is not only affected by the thick- 
ness of the body but is influenced essentially by the 
deviation of the stream density along the stream lines 
from its maximum value at the sonic velocity. 

For axisymmetric flow, the deviational ordinates of 
the stream lines tend to zero for large distances, but the 
area between the deformed stream line and the line y = 
constant at which this stream line starts also go to in- 
finity as y increases. 

The question can be raised as to whether the pattern 
obtained is of significance also for flows with free- 
stream velocities slightly different from sonic. It is 
probable that the flow picture obtained here can be ex- 
pected at points in a field which are not only sufficiently 
far from the body that details of its shape are unim- 
portant but which are also remote from portions of the 
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TABLE 1 


Table of Functions (See Figs. 3 and 4) 


= Axisymmetric 


+ ¥ 
t f f} = 
g 7 yo" 
94 —1_ 26124 —(). 24929 0.00572 
9 2 —1.31395 —(0). 27854 0.00704 
Ri - 1.37298 —(0.31257 0.00878 
at 3 - 1.439338 —(0.35186 0.01118 
_16 ~1.51408 —0. 39649 0.01452 
ha —1. 59824 —(). 44572 0.01928 
= 9 — 1.69255 —(). 49749 0.02607 
=) 0 —1.79717 —(). 54799 0.03567 
_0 8 -1.91131 —(0. 59168 0.04905 
~0.6 —2. (03294 —().62184 0.06730 
=().4 2. 15856 -().63118 0.09116 
=().2 2. 28347 —(0.61284 0.12123 
0) 2 40146 —0.47152 0.15716 
0.2 2.50539 —0. 505389 0.19914 
0.4 2.58729 —0. 34082 0. 24466 
0.6 2. 63856 -~(). 16483 0.29143 
0.8 2 65012 +). 05695 0.338596 
1.0 —? §1260 +0. 32653 0.37328 
1.2 2.51664 0.64506 0.39780 
1.4 2.35221 1.01584 0.403835 
1.6 2.10859 1.43541 0.38122 
1.8 1.77678 1.90599 0.32209 
2.0 — 1.34552 2.42768 0.215388 
99 —() 80487 3.00104 0.04936 
2.24 0.71213 3.11912 0.02000 
2.24 0.71213 0.15278 0.02371 
2.27 —(). 71038 0.15160 0.02078 
2.3663 0.69141 0.18900 0.01797 
9 554 —0. 66903 0. 12527 0.01468 
9 7735 0.64146 0. 10985 0.01149 
3.113 —(). 60477 0.09159 0.00811 
3.795 —(). 54724 0.06753 0.00447 


flow field which are strongly influenced by the changed 
free-stream Mach Number. To be more specific, for 
high subsonic velocities the sonic line will first go away 
from the body and then curve downstream to end in the 
shock that terminates the supersonic region. The pres- 
ent solutions would be applicable in a region that is not 
only sufficiently remote from the body but which also 
is under the influence of Mach waves that originate 
along that part of the sonic line which has a shape 
given by the present solution. 

In the supersonic case, the sonic line will end at the 
shock that starts in front of the body. The present 
solution holds for that part of the flow which is under 
the influence of Mach waves that start at a portion of the 
sonic line not too close to the shock and not too close to 
the body. Obviously the range of Mach Numbers for 
which the present solution gives an approximation is 
not too wide. 

In this connection, another fact may be mentioned. 
As can be seen from Fig. 2, the velocity behind the 
shock is close to sonic. At subsonic free-stream Mach 
Numbers, the transition from supersonic velocities to 
This 


shock does not lead back to the free-stream Mach Num- 


subsonic velocities occurs by means of a shock. 





— Planar 
Fe f f? y /5y 
—3.4 5.53117 —().78014 0.058211 
—3.2 5.37300 —(). 80023 0.06262 
—3.0 5.21069 0.82156 0.06741 
—2.8 5.04899 —() 84418 0.07268 
—2.6 t 87264 0.86810 0.07851 
—2.4 +. 69639 —0. 89326 0.08503 
—2.2 $.51499 —().91949 0.09239 
-2 (0 +. 32827 0.94647 0.10078 
—1.8 1.136138 0.97361 0.11042 
-1. 6 3.93864 0 0.12159 
—1.4 3.73609 l 0.13462 
‘2 3.52919 — 1.04298 0. 14992 
—1.0 3.31921 1.05402 0 16791 
O.8 3.10821 1.05236 0G. 18908 
—0.6 2. 89929 — 1.03217 0.21384 
—().4 2.69683 0.98655 (0). 24247 
=0.2 2. 50667 —(). 90793 0). 27499 
0 2.33618 0.78864 0.31098 
0.2 2.194238 0.62142 0. 34942 
0.4 2.09108 0.39994 0. 38859 
0.6 2.03813 0. 18859 0.42567 
O.8 2.04782 +(). 22693 0.45704 
1.0) 2.13333 0.64000 0.47787 
1.2 2.30832 1.12198 0.48179 
1.4 2 58735 1.67840 0.46227 
i.6 2 95804 2 30837 0.41059 
| .3 3.51628 3.01346 0.31714 
2.0 4.19618 3.79470 0. 17099 
2.03 4. 20818 3.80737 0.13804 
2.03 4. 20898 1. 36064 0.13058 
3.1666 5.50417 0.93312 0. O81L09 
1 749 6.80248 0.73246 0.041538 


ber, but it will lead to a Mach Number that is close to 
the free-stream Mach Number, probably even below 
(If the shock would 
lead to a Mach Number above the free-stream Mach 


the free-stream Mach Number. 


Number, then the stream lines behind the shock must 
diverge along the axis of symmetry, since the stream 
density under free-stream conditions would be lower 
This 


From this idea, the approxi- 


than the stream density right after the shock. 
appears improbable.) 
mate location of the other shock can be found, and it 
can be seen that, even at rather high subsonic Mach 
Numbers, this shock will lie relatively close behind the 
body. 

If one considers instead of a closed body a half-body 
with a finite thickness at infinity, then the asymptotic 
flow field will not be changed essentially. (In contrast 
to incompressible flow, the effect of the finite thickness 
at infinity could probably be investigated by means of a 
perturbation theory.) In this case, there is a marked 
overexpansion out to distances far from the body—.e., 
in spite of the fact that the contour of the body in its 
upper half has at no point a negative slope, the stream 
lines in the flow field will have points with a negative 


slope. 
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A Remark on the Stability of the Laminar 
Boundary Layer in a Compressible Fluid+ 


Martin Lessen 

Professor of Aeronautical Engineering, The Pennsylvania State 
College, State College, Pa. 

March 30, 1953 


iy A RECENT PAPER, Lin! discussed the stability of the super 
sonic boundary layer to small disturbances with the conclusion 
that, under certain conditions, a three-dimensional disturbance 
could be more destabilizing than a two-dimensional disturbance. 
It follows that, for such a situation, one must, in effect, investi- 
gate the stability problem for a range of three-dimensional dis- 
turbances in order to find the most destabilizing one. 

It seems that it might be possible to consider such two-dimen- 
sional disturbances that the minimum Reynolds Number for neu- 
tral stability thus obtained would fall in an actual region of cer- 
tain stability. 

With a fixed dimensionless profile, consider the minimum neu- 
trally stable Reynolds Numbers defined as follows: 


R = F(M, B) 

R° = F(M,0) 

R = F(M, 0) 
where .)/ is the free-stream Mach Number, 8 is the angle between 
the direction of the disturbance propagation and the principal 
flow direction, R = R cos 6 as given in Squire’s transformation,? 
and \f = M cos @ as given in the generalization of Squire’s trans- 
formation by Dunn and Lin.* 

At sufficiently high Mach Numbers, Lin! found that there ex 

ists a B* such that R* = F( AM, B*) < F(.M,0) or R* < R° 


Now 
R* > R* = R* cos B* 
R* = F(M*, 0) 
where \/* = V/ cos B*, or R* < R° and M* < M for the same 


dimensionless profile. At 1/ = 0, by the Squire transformation, 


R = ReosB = R® 


It therefore seems that, at large values of ./, R R°® and, at 
small values of 1/7, R = R°. 

Keeping the same dimensionless profile, it therefore seems, by 
continuation, that a minimum R will be obtained for 6 = 0 and 
M=0 

It therefore seems, although not rigorously, that certain stability 
to small disturbances can be obtained using the last equation of 


that is, if the disturbance is incompressible. 


reference 4—namely, 

U U - 
~ ¢) ((Bo')’ — aXpv)] — (BU')'v = — i 
( U pi a*( pr p U aR. 


t This work carried out under ONR Contract NONR-656(01) 


500 


where R, is defined by the viscosity at 1’ = c and the balance of 


the notation is given in reference 5 
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Notes on the Phugoid Damping 


Frank W. Heilenday 
Cornell Aeronautical Laboratory, Inc., Buffalo, N.Y 
April 1, 1953 


—— USUAL AIRCRAFT PHUGOID OSCILLATION consists of an 
altitude, air-speed, and attitude variation of almost neutral 
stability. In order to maintain a stable flight path during exact 
ing maneuvers, the pilot (or autopilot) must then react to damp 
out the phugoid motion caused by air gustiness or elevator mo 
tion. Since the phugoid period can be estimated to be one-fifth 
of the true air speed in miles per hour, this long-period oscillation 
can normally be controlled easily by the pilot. Whenever the 
pilot cannot devote constant attention to his flight attitude (as 
in the case during instrument flight or gunnery runs), the low 
damping of the phugoid may become troublesome, leading to 
pilot fatigue and to aiming error. In the case of a modern high 
speed fighter, this oscillation becomes of such a long period as 
to appear as an insidious divergence in attitude and air speed 

Fundamental research into the importance of the phugoid in a 
pilot’s evaluation of an aircraft is now being performed at the 
Flight Research Department of the Cornell Aeronautical Labo 
ratory, Inc. Two aircraft, an F-94 and a B-26, have been instru 
mented to allow a wide range of phugoid dampings and periods 
to be test-flown. A small auxiliary pitching surface, installed in 
each airplane, will be positioned by signals from air-speed changes 
and rate of change of air speed to provide C,,,, and C,,; artificial 
derivatives. The effect of these added derivatives can be seen 
in the phugoid phase diagram in Fig. 1 for the F-94 airplane 
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It can be seen by the increased phase lead between “@ and u 


that the derivative C,,;, added 50 per cent critical damping to the 
phugoid. The phasing of @, the aircraft attitude, relative to 
V, the nondimensional air-speed increment, has been decreased 


ul 
This change in 


from the normal phugoid oscillation, however 
the nature of the phugoid must be considered in the evaluation, 
is well as the desired increased damping 

The use of C,,,, to decrease the period of the phugoid is seen to 
increase the magnitude of the attitude change relative to the air 
speed change as shown by the increased length of the @ or atti 
tude vector. The combined variation of both C,,,, and C,,; in 
flight will be used to evaluate many different combinations of 
optimum desired and minimum tolerated period and dampings 
A variation of period from 50 to 150 sec. and of damping from 70 
per cent critical unstable to 70 per cent critical stable can be tested 
on flights of the F-94. 

Although no flight results have as yet been published, the re 
ports in references 1 and 2 have been issued on the theoretical 
investigations. These studies are being performed for the Aero 
nautical Research Laboratory of WADC under the guidance of 


P. P. Cerussi 


REFERENCES 
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On the Nonexistence of Four Confluent Shock 
Waves 


Julian L. Davis 
Republic Aviation Corporation, Farmingdale, L./., N.Y. 


March 26, 1953 


I STUDYING THE STATE OF AFFAIRS existing locally about the 
point of intersection of four shock waves, it was discovered 


FORUM 501 


that it is impossible for these waves to coexist under certain con- 
ditions. Specifically, consider a supersonic flow field in which 
there are four shocks (.S;, S», S;, S;) confluent at ? which divide 


the field into the regions shown in Fig. 1. Let the shock con- 


figuration be stationary and let the flow move from left to right 
so that region O is undisturbed, () and (©) are once disturbed, and 
3) is twice disturbed. Suppose the flow in each region is uniform 
(region 7 has pressure pi, density p;, etc.) and the four shock 
strengths, as measured by their pressure ratios, are unequal. It 
is further assumed that each shock may be represented by a 
straight line of zero thickness. It shall be shown that this situ 
ation is impossible—i.e., four confluent shocks of unequal strength 
cannot coexist if the flow in each region is uniform 

The proof given here is based on the one applied to the steady 
three-shock or Mach configuration.! The starting point is the 
set of Hugoniot expressions which shall be used in the form of 
shock strength (pressure ratio) as a function of density ratio for 
any two adjacent regions 7, 7. 


Wiy = (Vij w)/(1 Brij (1 
= pi/pj, rij = pi/pj are the pressure and density ratios 


The 


where 7j; 


ind yw = (¥ 1)/(y + 1), y being the ratio of specific heats 


identity 32713 = mo702 then becomes 


Geers Cer (or Ce 


(2). 


A given gas fixes uw (say w), and so w must be a root of Eq 
, and for a physically real 


But uw = 0, +1 are also roots of Eq. (2 
izable gas up ~ 0, +1 Furthermore, it is seen that Eq. (2) isa 
This leads to a contradiction, meaning that the 


cubic in yu 
The proof has 


physical setup leading to Eq. (2) is impossible 
thus been established 
Suppose shock .S; equals S: 


G mi & mh ; 


Then 327; = 1 or 


which becomes (on using ’32%13 = "i2) 
(1 rio)? = | r (3) 
But riz = 1, since by the perfect gas law ra) = my /7 2 and 
Tx = 72/T, = xn = 1 
because S; = Sv. Eq. (3) is an identity and, hence, holds for 


B= 4B; In fact, the case S; = 
Ss is equivalent to a classical reflection problem.” 

It might well be that four unequal confluent shocks can exist 
3) from P Cour 


there is no apparent contradiction 


with a line contact discontinuity D in region 
ant-Friedrichs! postulates a contact discontinuity line for the 


three-shock configuration 











@ 





Fic. 1. Four-shock configuration 
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Evaluation of the Momentum Integral 
Equation for Turbulent Boundary Layers; 
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F’ THE PAST FEW YEARS, workers in the field of turbulent 
boundary layers have recognized that the values of wall shear 
stress calculated from the von Karman integral momentum 
equation may be seriously in error, especially in adverse pressure 
gradients and near separation. Newman! and Wallis? in Aus- 
tralia and Bidwell? and Rubert and Persh‘ at the NACA have 
shown that an appreciable correction to the mean-flow equation 
results from retention of the turbulence normal stresses in the 
boundary-layer equations. Goldschmied® has integrated the 
Reynolds equations without making any of the boundary-layer 
assumptions and has shown that the turbulence contributions to 
the momentum thickness should not be neglected. In a recent 
note, Van Le® has also derived an integral momentum equation 
that includes the turbulent fluctuation terms. Unfortunately, 
his derivation, as published, contains an error that affects the 
order of magnitude of one of the turbulence terms. The writer 
has also carried out an integration of the complete Reynolds 
equations,’ obtaining a result equivalent to that of Goldschmied® 
but retaining the usual definition of the momentum thickness 
and expressing the deviations from the von Karman momentum 
equation as four correction terms. In the report,’ the correction 
terms are evaluated as to order of magnitude, and it is shown that, 
except at separation, only the correction term involving the 
turbulence normal stresses is important. The resultant approxi- 
mate equation is 


dé T, 260 + 5* dp I fc o 


/ att ; 1) 
dx pu;? pu? dx uy? J9 Ox (u F ) dy 


all higher order corrections being omitted. 

By relating the anisotropy of the turbulence to the shear, an 
approximate expression for the correction term in terms of mean- 
flow quantities can be evaluated, as follows: 

u’? — y’2 = —2u’v'/tan 2a (2) 
where a is the angle which the principal axis of the turbulent 
stress tensor makes with the direction of the mean flow. Taking 
the derivative outside of the integral, the correction becomes 


> € 6 n'v' 
. - — pu’ ‘ 
CG = f = & (3) 
pu? dx J9 tan Qa ~ 


Qa 


Dryden’ has suggested that the turbulent shearing stress may be 
related to the local nonisotropic turbulence and that a may be 
approximately constant. Examination of the data of Schubauer 
and Klebanoff® and of Laufer" tends to confirm this hypothesis, 
so that, at least for the purposes of evaluating this correction 
term, we may assume tan 2a to be constant with the numerical 
value ~0.75. Evaluation of the correction term then reduces 
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to calculation of the integral of the shear. Granville'! has showy 


that this integral is approximately proportional to the displace. 
ment thickness, so that, finally, 


: 27 d 6 2.7 d an 
C; = ; (—pu'v’) dy = : J, (0:0068* 1?) (4 


pu,? dx u,? d: 


Expressing 6* as 06H and taking the derivative with respect to x 
yields the corrected integral momentum equation in the form 


Te 6 dp 


tf 
+ (2 + 0.97H) + 0.0160 


dé pu,” pu? dx dx (5 


dx 1 — 0.016H 


The importance of the turbulence correction is realized when one 
notes that the factor in the denominator gives a 2 per cent cor 
rection even for flows with zero pressure gradient. The correc 
tion can amount to about 10 per cent of d6/dx for flows in strong 
adverse pressure gradients. 

Approximate evaluation of the other correction terms in similar 
form has been carried far enough to show that, except possibly at 
separation, the errors involved in the experimental evaluation of 
the terms in the momentum equation and the errors in the assump 
tions used to evaluate the major correction overshadow these 
other terms. It is interesting to note, however, that use of the 
static pressure gradient measured with wall piezometers, instead 
of the free-stream pressure gradient, reduces these higher order 
corrections appreciably. 
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yg AN INTERESTING PAPER on hypersonic viscous flow over 4 
flat plate,’ Lees mentioned that for the region of small longi- 
tudinal distance x from the leading edge of a flat plate, the Pohl- 
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hausen parameter A is approximately constant. In this note, 
the author discusses the conditions under which A(x) = constant 
and also the boundary condition of pressure at the outer edge of 
the boundary layer in hypersonic flow over a flat plate. 

If the pressure p at the outer edge of the boundary layer may 
be expressed as 


p = p(di/dx)’ (1) 


where 6 is the thickness of the boundary layer and pp) and are 
constants, then 


A (x) = constant = A (2) 


for an insulated plate when the Prandtl Number is assumed to be 
unity. 

The proof is as follows: 

From the boundary condition at the surface of the plate if 


u/s, = f(y/5) (3) 
we have 
6° dp 6? dp 
-A = = (4) 
M,, U2 dx MM, Uy dx 


where y,, is the coefficient of viscosity at the wall, uw. is the ve- 


locity at y = 6, and ™ is a constant, say #2, atx = 
Now, if A (x) = Ao, Eq. (4) gives 


, (di\"” 1 q25 By Uy A rp 
52 - = — = constant (9) 
dx dx? np 
Hence, 
é a itt) (n+2 (6) 


where 6p is the proportional constant, and then 


o n-+1\' 
P = § g /F TS) « ) (7)* 
po n+2 


Substituting Eqs. (2), (6), and (7) into the momentum integral 


) x 
f 62 = F(A) f, G(A) dx (8) 
Pr 0 


where F(A) and G(A) are given functions of A, the resultant equa- 
tion is independent of x; thus, A(x) = Ag is satisfied in the pres- 
Eq. (8) may be used to determine the constant 


equation t—i.e., 


ent problem. 
\ 

The results of Lees and others” 
author will give examples for the 


> correspond to the case m = 2 


In the following section the 
cases other than » = 2 and also discuss the cases when A (x) is 
not equal to a constant 
The boundary condition of pressure at the outer edge of the 
boundary layer is obtained from the oblique shock relation 
p 27 M_,? sin? 0 y-1 
= _ (9) 
p.. y+ 1 ¥+1 
Where 6 is the shock angle with respect to the uniform flow of 
Mach Number / in front of the shock. In order to satisfy 
the asymptotic condition at x = ©, we write in the problem that 
a flat plate is placed in the direction of a uniform plate (Lees’ 
problem of reference 1) 


@=a+t+, (10) 


where a = sin~!(1/M_,)is the Mach angle of the uniform flow 
M and 6; is the shock angle induced by the boundary layer 
flow. If a + % is small, Eq. (9) becomes 

* It was pointed out to the author by Dr. S. F. Shen that the present con 
dition is similar to the condition for the existence of similar solution in the 
theory of the boundary layer of an incompressible fluid—i.e., if pressure is 
Proportional to a power of x, the solution is similar and A(x) is a constant.* 


Tt See Eq. (19) of reference 2 
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Fie./, Hypersonic Viscous Flow Over A Wedge At An 
Angle Of Attack 


p 4yM_.41 27 M?2_ 4:7 
ei + + (11) 
P-« ¥ +f y¥+1 
The value of 0; may be determined by the transverse velocity com- 
ponent v2 at the outer edge of the boundary layer and the oblique 


shock relation 


v, _ M?_,, sin 20 — 2 cot 0 (18 
= e) 
Ue 2 + M?_..(7 + cos @) 
The result is 
v2 - dé 
= F(a):0; = (13) 
Ue dx 


where #» is the longitudinal velocity component at the outer edge 


of the boundary layer. F(a) is a given function of a. Only 


when 6; > > a 
F(a) = 2/(7 + 1) (14) 


Combining Eqs. (11) and (14), when @; >> a, we have Lees’ ap 
proximate boundary condition of reference 1 
For those ranges of x, where the term J/ 
gible, the boundary condition (11), at least the last two terms on 
should be used with relation (13). In this 


, 9; is no longer negli 


the left-hand side, 
case, there will be no similar solution, so that A(x) is not a con- 
stant. The exact solution may be obtained by the conventional 
method of series expansion 


Another interesting result from Eq. (11) is that, if 
M_. 6 << 1, F(a) = 4/(y7 + 1) and 
p/Po = 1+ 7M, (di/dx) (15) 


which agrees with Lees and Probstein’s first-order boundary 
condition for large longitudinal distance,® even though their 
considerations are different entirely from the present analysis 

If we consider the hypersonic flow over an insulated flat plate 
or wedge at a large angle of attack (Fig. 1),{ on the expansion 


side, we have approximately 


-1 15 \27/(9-1 
- -(7 — My. ; ) (16) 
> +20 - ax 


Eq. (16) satisfies the condition of Eq. (1), so that the similar 


solution exists, A(x) = constant and 


6 ~ x87-D/4y—2) (17) 


8 
for y = 14,56~x’"° 

t For convenience, the expansion and the compression side may be treated 
separately with different reference axes as shown in Fig. | Ifao=A 0, 
case of reference 1, which is shown in the 


then 4% a and we have Lees’ 


compression side 
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On the compression side, we may write 


6=h%+h (18) 
where @ is the shock angle without the induced effect and 
p p 2y71\/?_. 0: sin 26, 
p p 1 ine a ; 
24 MJ? 6,;?_.. cos? % 140) 
y+ 1 


Here the last term is generally negligible in comparison with the 
other two if 6) is large and 1/7... is still much larger than one 
Only when the second term is much larger than the first term, 
similar solution, A(x) = constant, may exist, and then 


6~mx/3 (20 
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SUMMARY 


The effectiveness, as a control surface, of a variable incidence slender wing 
on a cylindrical body, is determined The lift of such a configuration is com 
pared with the lift of a wing-body combination having both components 


at the same angle of attack and with the lift of an isolated wing 


SYMBOLS (PARTIAL LIsrT) 


v, 9,2 coordinates in streamwise, spanwise, and vertical directions, 
respectively (Fig. 1) 
a body radius 
spanwise distance from origin to maximum wing span 


c complex variable (y + iz) 
q dynamic pressure (1/2 pl?) 
Lat lift of wing-body combination (excluding forebody lift) when 


wing is at angle of attack a and body is at zero angle of attack 
Laa lift of wing-body combination (excluding forebody lift) when 
wing and body are both at angle of attack a 
La lift, at angle of attack a, of isolated wing formed by joining ex 
posed wing panels of given wing-body combination 


DISCUSSION 


pee PROPOSED MISSILE DESIGNS have incorporated vari- 
able incidence wings that serve both as lifting surfaces and as 
control surfaces. The wing panels are hinged rather than solidly 


connected to the central body. The flow about such a configu 
ration may be expressed as the linear superposition of the follow 
ing two flows: (a) the wing and body at the same angle of attack 
and (b) the wing at angle of attack and the body at zero angle 
of attack. 

The solution of (a) has been presented by Ward! and others, 


subject to the limitations of slender body theory. The solution 
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of (b) defines the effectiveness of the wing as a control surfacg 
and is the subject of the present note. Slender body theory js 
assumed toapply. The effect of the gap that may exist betweey 
the wing and the body is neglected herein and will be treated jy 
another paper. 

Consider the body to be cylindrical and at zero angle of attack, 
while the wing is at angle of attack a to the main stream L’ (Fig 
1). According to slender body theory, the cross flow in each 
ys plane may be considered independently. The boundary condj 
tions in the ys plane are that the normal velocity at the body sur 
face is zero, the normal velocity at the wing is w = al’, and the 


complex velocity v — iw behaves like 1/¢? for large ¢ The solu 


tion for the cross flow may then be shown to equal 


al ( “) 
9 I -9 . 


(w+al) = 


tcot | (a/y: 


In 
( Zaye ) ( “y( =) 
l 4 < yo + y + ] 
V2? a? \ : ye . y 
$ cot | (a/ye 
- ) 
"ah ig a? 
ve + y+ l 
\ y y 
For a point on the body (¢ ae’ 
q al 2A 
eo =a) == — (1 3S om 
( Zaye ) ( =} . 
- Yo + (2a cos 6) 2 ] l 
Vo? a? \ 3 Ve 
In 


( 2a Vo ( =) . 
‘ yo + (2a cos 6) 2 l 
2? a? \ ‘i y 


$cot '(a V2) 


a*\? 
Yo 4 (2a cos 6) ? l 
\ ys 


Following Ward,' the lift of the wing body combination may be 
determined from the asymptotic form of Eq. (1) as ¢ becomes 
infinite, and equals 


Lat 


2rqa 


2 a? a a \? 2a a? 
1 + cot ! 2 l } 

T y2” V2 y2 T V2 ye") 
Eq. (4) defines the wing effectiveness. It may be compared to 
the lift on an isolated wing, at angle of attack a, formed fron 
the exposed wing panels of the wing-body combination. The 


lift on such a wing is 
La/2xqa = y-?(1 (a/y2) |? 5 


Eq. (5) represents the lift that the wing-body combination would 
have if the body acted as a reflecting surface and carried no lift 
The ratio Lao/La is plotted as a function of a/ye in Fig. 2. Its 


value increases almost linearly from 1 at a/y. = 0 to 2 at a/y¥2 = 
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Slender wing at angle of attack @ in respect to cylindri 
cal body 
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Fic. 2. Lift effectiveness 


1. Thus, for the case of small wings on a relatively large body 


a/v. = 1), the wing has twice the effectiveness that would be 
computed by considering the body to act merely as a reflecting 
surface without carrying lift. Inspection of Eqs. (2) and (3) 


indicates that, for a/y. = 1, half the lift is on the wing and half is 
on the body 

Eq. (4) may also be compared to the lift that would exist if both 
the wing and body were at angle of attack a. The lift in this case 
is! 


Laa/2rqa = y2* [I (a?/y2?) |? (6 


) 


The ratio Lav /Laa is also plotted in Fig. 2 and has the value 1/2 


ata/y. = | Ward, in his discussion of Eq. (6), pointed out that, 
for a/y: & 1, the lift on the wing corresponds to twice the geo 
metric angle of attack of the wing-body combination (due to up 
Wash induced by the body) and that there is an equal amount of 
lifton the body 


to the case of any arbitrary orientation of the wing in respect to the 


The present investigation generalizes this result 


body. Thus, for a/y. = 1, the lift on the wing corresponds to an 


angle of attack equal to the sum of the geometric angle of attack 
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of the wing plus that of the body, with an equal amount of lift 
appearing on the body 
It might also be mentioned that Lay is equal to the geometric 


mean of Le and Lea to within 4 per cent That is, 


0.96 « La V Lelea < l (7 


Lagerstrom and Graham? have discussed some of the limita 
tions of slender body theory. In particular, they point out that, 
if there is no afterbody and the flight Mach Number is suffi 
ciently above 1, slender body theory will tend to overestimate the 
lift carried by the body. In this case it may be more accurate 
to consider the body as a reflecting surface. Such considerations, 


however, are beyond the scope of the present note 
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A Remark on ‘‘On Asymptotic Solutions for the 
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Sweden 


T° AN INTERESTING CONTRIBUTION to the Readers’ Forum, * 
Schuh deduces a differential equation for the asymptotic 
the the 


boundary layer 


solution of dimensionless temperature @ in thermal 


do £2 do 


é 
dé? 2 dé 


with € as a parameter, € > 0. The boundary conditions are 


E= QO, 6=1 
E> ow, a@— 0 


Schuh remarks that he did not find a solution to this equation in 
closed form and that a numerical solution has not been attempted 
by him. 

Perhaps, it may be of some interest to point out that to con 
struct the required solution by using purely numerical methods 
would be a difficult task because the general solution of the equa 
tion is highly singular at infinity. Fortunately, this is not neces 
sary. In fact, the required solution may be expressed in closed 


form 


oe ee) ler (-a) 


The function I'(#) is the gamma function and the function ,/i(a, 
c, 2) is the confluent hypergeometric function, defined as a solu 


tion of Kummer’s differential equaticn 


dy di 
t+ (¢ ) av 0 
dz? d 
* Schuh, H., Readers’ Farum, JOURNAL OF THE AERONAUTICAL SCIENCES 


Vol. 20, No. 2, pp. 146-147, February, 1953 








506 JOURNAL OF THE 


A Positive Method for Locating the Zero-Order 
Fringe in Interferograms Produced with 
Filtered Light 


Henry W. Wedaa and Edward W. Price 
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I* EXPERIMENTAL RESEARCH utilizing the Mach-Zehnder inter- 
ferometer, it has been extremely difficult to trace the zero- 
order fringe in filtered spark interferograms by visual means. 
This is a serious disadvantage, since a particular fringe cannot 
be identified in a sequence of pictures as is ordinarily required in 
calibration of the optical system.! One method of locating the 
zero-order fringe by the use of a photometer has been discussed 
previously.” 

In the present method, using a magnesium spark and a suit- 
ably placed interference filter, the zero-order fringe is made easily 
identifiable provided the interferometer is so adjusted that the 
color dispersion in each branch of the interferometer is approxi- 
mately the same. The interference filter is mounted in front of 
the objective lens of the camera at a right angle to, but still in the 
parallel beam of, the light coming from the interferometer (Fig. 1). 
It is mounted so that most of the light passes through the filter 
and a small portion does not. The resulting interferogram (Fig. 
2) will have one part containing many interference fringes due to 
the filtered light and one small strip containing only a few fringes 
due to white light. If the interferometer has been properly ad- 
justed for white-light fringes * and the dispersive media in the 
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Fic. 1. Illustration of position of interference filter with respect 


to interferometer and camera. 
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Fic. 2. 


Interferogram with part of field unfiltered to indentify 


zero-order fringe 


branches of the interferometer are matched, the zero-order fring 
can be traced from the unfiltered part of the interferogram (wher 
the contrast between the zero-order fringe and the remaining 
fringes is great) through the filtered part of the interferograt 
(where the contrast between the zero-order fringe and the x 
maining fringes is such that the zero-order fringe cannot be identi 
fied by the eye). Thus the zero-order fringe in the unfiltere 
part of the picture is used to label the zero-order fringe in the 
filtered part of the picture, and this is accomplished without auxil 
iary apparatus or operations. In practice, all that is required is 
that the interference filter be provided with a narrow band that 
is uncoated, so that a narrow strip in the interferogram is phot 


graphed with white light. 
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The Relation Between Open-Loop Stability 
and Frequency Response Diagrams for Linear 
Systems* 


Daniel O. Dommasch 
Chief Engineer, Test Pilot Training Division, NATC, U.S. Naval Air 

Station, Patuxent River, Md 
March 20, 1953 


(1) INTRODUCTION 


6 iar PROBLEM OF DETERMINING whether a mechanical, electri 
cal, hydraulic, or combination system describable by linear 

* The mathematical developments presented herein were suggested an¢ 
checked by Daniel C. Lewis, Professor of Applied Mathematics, The Johns 


Hopkins University 
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differential equations is stable resolves itself into determining 
the nature of the roots of its characteristic equation. The char- 
acteristic equation may be obtained either by assuming that the 
solutions to the differential equation set are expressible in the 
form of a linear combination of exponentials Ae*t (where the A’s 
are constants determined by the initial conditions of the problem ) 
or by using operational methods. Whatever method is used pro 


yides an algebraic equation of the type 
. + ° + k, 12 
es 


s* + kys"—! + kes" + kn = 0 (1) 


where the coefficients k; . . are real and the roots in general 


are complex. 
If the characteristic equation has no roots with positive real 
parts, the system is stable; otherwise, it is unstable. 
It is the purpose of this paper to demonstrate that a plot of the 
quantity w = 1/characteristic polynomial for pure imaginary 
to —7a 


yalues of z from +7 may be used to directly deter- 


mine if an open-loop system is stable. For many systems the 
plot, as described above, may be shown to be the same as the 
system’s frequency response diagram, in which case the function 
w defines the system’s ‘‘vector gain’’? or ratio of steady-state out 
put to input for a sinusoidal forcing function of frequency y. In 
this paper a plot of the function w for any system is referred to, 
for want of a better name, as the frequency response diagram, 
even though this usage is not precise for all systems 

The procedure is similar to the Nyquist! criterion presently used 
to determine the characteristics of closed-loop systems and, like 
the Nyquist and other stability criteria, is based on the Cauchy 


equation 


1 fd log w 
| og — 
Qri dz 


(sum of orders of zeros 2 


sum of orders of poles) (2 


where w = f(s) is an analytic function of the complex variable z 
and the zeros and the poles referred to in the formula are those 
zeros and poles of w which are surrounded by c, a closed contour 


inthezsplane. It is assumed that w has no zeros or poles on ¢ 


2) THE SEMIGRAPHICAL OPEN-LOOP STABILITY CRITERION 
To utilize Eq. (2), it must first be transformed slightly. This is 
wccomplished by noting that 
log w = log | w)| + 164 
where w is the modulus of w and Oy is its argument. Thus, 
| °*d log w (Oy 
dz = (log |w| + vw). = (3 
2rie dz 2ri 2 


where (Oy). is the change in @y as s makes one continuous circuit 
of the boundary c in the positive sense. 


Combining Eqs. (2) and (3) 


wv )-/2% = (sum of orders of zeros — sum of orders of poles 


(4) 


To utilize Eq. (4), consider that f(z) is the left-hand side of the 
characteristic equation—i.e., 
Ros" 2 +. t kris +k, (5 


is) =" + 


b ont —1 
15 T 


Since f(z) is a polynomial, it has no poles although it has m zeros 


orroots. Consider next that 


f(z) (6) 


w= 1 


since f(s) has no poles, w its reciprocal can have no zeros but only 
poles corresponding to the zeros of f(z). Also if Rk, is not zero (it 
is never zero for oscillating systems), then w is finite for z = 0 
Now, the problem is to determine whether f(z) has any zeros 
roots) in the positive right half-plane which, in view of Eq. (6), 
amounts to determining if w has any poles in this half-plane 
The problem can be answered by mapping the entire right half- 


plane contour onto the w plane and counting the number of 
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revolutions the vector w makes as z follows the boundaries of the 
right half-plane contour. This is similar to the Nyquist proce 
dure for the closed-loop system; however, there is an important 
difference—namely, all the circuits nade by w will not be discern- 
ible directly as in the Nyquist case. To see why this is so, con 
sider Fig. 1. 

Allow z to go from 7R to —7R and thence around the border of 
the large circle shown in Fig. 1. 
©, so that at this limit, w —~ 0, and there 


Examination of Eq. (5) reveals 
that f(s) ~ ~asR—> 
fore the entire semicircle contour is mapped onto a small neigh- 
borhood of the origin of the w plane by the transformation 
Thus, speaking somewhat roughly, mapping the boundary of the 
infinite half-circle of the s plane into the w plane through the 
transformation (6) entails evaluation of w only for imaginary 
+1 Now, al- 
though points in the s plane where |s| = © are all transformed 
to the origin of the w plane, there still remains a rotation at the 
visual examination of 


values of s extending from z = tos = i 


origin which cannot be detected from a 
the plot in the w plane (obtained from the transformation of the 
infinite half-circle of the z plane). The reason for this is that ro 
tation of the vector w can still exist even though the modulus of 
the vector approaches zero. It follows that, to apply the Cauchy 
equation [Eq. (2)] to our problem, we must first evaluate the in 
visible rotation at the origin of the w plane. This evaluation is 
accomplished as follows 


Referring to Fig. 1, the problem is to determine 


d log % 
1(R) - {| d: 
(R dz 
as 0 goes from —7/2 to 2/2. Now, by definition, w = 1/f(s 


therefore, log w = —log f(s) and 
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i d log f(z) P 1 d{[f(s)| 
I(R) = - dz = — : dz and 
c(R) dz c(R) f(z) dz 
Since 4 
f(s) =o" +h" '+...4+ 8, sill | 
P P no"-1 + k(n —1)s"-2?+. + Ry 4 F Y, f 
I(R) = - dz me 
c(R) "+k '+...4+h, 
‘ TI 
factoring n/z from the numerator, dividing numerator and de- fF , eque 
. . ° 
nominator by s", and letting { oe port 
| b Path of w for 2= tinteniar \ that 
14 n k , n—I a ‘ Note the two complete counter 
t ne ee ee oe no"-! ie clockwise rotations of « about is al 
F i \ the origin e 
(s) = g0° ‘ - \ x ence 
1 + (ky/z) +... + (Rn/s") ¥ ) w 36.25.) ~ k 
a ww O° right 
[(R) becomes ie ! 9 
t| - he 
. a FIG. 28 the 
I(R) = f - F(z) dz STABLE QUARTIC curv 
c(R) ¢ 
. \ , not 
Now, z = Re’ and dz = iRe*’ do + e'’ dR, where R = mod. z and J \ , real 
@ = arg. cs. Moreover, if we integrate only over the contour i 
“X‘ 
of the semicircle dR = 0, so that . 
————— 
w/2 9 90° 4 “¥ 270° In 
as W( Ret @ is the vector from the origin 3 
I( R ) t — 7/2 n I ( Re } de to Gny point on the frequency servi 
response curve resp 
Therefore, noting that from the definition of F(s), 7 Note there are no rotations of = craft 
‘ weabout the origin _" 
. . . * 4 ‘ - 
lim F(z) = lim F(Re™) = 1 ils Geen Zs —" ‘\ tem 
Is|—< R— 180° | =——- od Sn ancgall exist 
i itl SSS = a © root: 
uniformly with respect to 6, 
= —— r 
°,/2 — 
x /2 270° , 
‘ : . } syst 
lim (J(R)) = -1 f n| lim F (Re'”’) | de = FIG. 28 yst 
-/2 FIG 28 
Ro w/e Ro UNSTABLE QUARTIC and 
. , ? 
- proc 
if n dé = ni (7) yee F : ; 
n/2 dividing into real and imaginary parts gives 10W' 
. os . eT ee ee ak ca) soa hae 7 : : ; ; for « 
Eq. (7) defines the “‘invisible’’ rotation at the origin fiiy) = (y* — 18y? + 16) — i(6y' D4y aie 
. . . o. . . . . ba = ~ - as > ( - 
Dividing the contour of integration in Eq. (2) into two parts, st ac : = : 
. ; : . i e6 laking the reciprocal of f(7y) and multiplying numerator and the | 
the first as evaluated above extending over the semicircle and the é : : Ks Seni . 
; : . é 7 denominator by the conjugate of the denominator, w is found to is al 
second extending from +/© to —7~™, we find that : : ig 
be origi 
l io d log w 
— . : ? smor 
dz mi} = Ce Sy? + 16 + i(6y 24y 
2ri tic dz w= : 
; : , (y' I8y? + 16)? + (6y 24y)? 
(sum of orders of zeros sum of orders of poles) (8) 
, — whence i) 
and utilizing this in Eq. (4) results in ary, | 
Z ia w| = mod.w = ; 
(Ow )s= +3 « n _ V (y' 18x? + 16)? + (6y 24y the 
2r 2 
by? 24y 
(sum of orders of zeros — sum of orders of poles) (9) 6 = arg. w = tan”! ‘ = : 
yt — 18y? + 16 
From Eq. (9), we immediately deduce the following rule: a sa , ; ; 
" : : F ; : rhe brief numerical calculations required to determine w, ane 
If a system is stable, mapping the imaginary axis of the z plane . . . ; 
nm tt lirecti f +i ‘ ; ; th f : @ are omitted here for brevity, but the resulting frequency re 
in the direction from 1© to —1© using the transformation, : . ws , 
8 sponse curve is shown in Fig. (2A). We note that, because the Th 


w = 1/characteristic polynomial, must provide a counterclock . , -— : : 
portion of the curve below the real axis is the mirror image of the 


al a 


wise rotation of w in revolutions about the origin exactly equal : ar f 
a : : portion above, it is only necessary to compute w for values of ) 
to one-half the degree of the characteristic polynomial. If one 


é tae i : ranging from ©+ to(Q. From an examination of the frequency 
considers the variation of w as z goes from —i@ to +i, the ro sie : KY 
‘ ° ae response curve, It 1s seen that, as s goes from +7:1© to 1a, . 
tation of w is reversed. eA ‘ : ; . . A 
a , a so . a w rotates about origin two times in a counterclockwise direction, er 
ro use this criterion, the characteristic equation is first ob- : cee ais rhe : - : C 
. a ; ; which satisfies the stability criterion described previously in the 
tained, and a map of the /y axis is drawn using relation (6). Mar 


pis : : ; paper. (Note that the rotation of w is given by the change 
The rotation of w as s goes from +7 to —/@ isthen notedto  . me 
: 5 fi in 0w as z moves from one limit to the other 
determine compliance with the above rule < : : 

Consider now an unstable system described by the character- 


istic quartic | 


(3) ILLUSTRATIVE EXAMPLES 





st +2 27 +4=0 root 
Consider the quartic characteristic equation : side; 
with roots 
f(z) = 24 + 623 + 18s? + 245 + 16 = 0 (-1).(-2), 0+0.01 cons 
ae bico 
which has the roots ; , ; ; : 
‘ Letting s = zy, this equation gives latic 
(—1 +8), (— 1 — 5), (—2 + 28), (— 2 — 23 (yi + +4) + Hy . 
w= —— : lar 
¥ . ae « 
and therefore represents a stable system. Replacing z by iy and (yt + + 424+ (y 
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and 


4)? + (y3 2y)? 


The frequency response curve determined by the foregoing 
equations is shown in Fig. (2B), and it is noted that the visible 
portion of the curve provides no rotation about the origin and 
that therefore the system is unstable. Referring to Eq. (9), it 
is also noted that the frequency response curve reveals the exist- 
ence of two poles of w [i.e., zeros or roots of f(z)] in the positive 
right half-plane. This follows from substituting » = 4 into Eq. 
(9) and recalling that w has no zeros by virtue of the nature of 
the characteristic polynomial. Thus, the frequency response 
curve serves not only to determine whether a system is stable or 
not but also to determine the number of roots having positive 


real parts. 
(4) CONCLUDING REMARKS 


In designing modern aircraft employing hydraulic-electrical 
servo systems, it is necessary to determine air-frame frequency 
response characteristics. It has been shown here that the air 
craft frequency response diagram (or such a diagram for any sys 
tem) immediately indicates whether or not dynamic stability 
exists. Moreover, for an unstable system the number of unstable 
roots is also given directly by this diagram. 

The procedure developed herein is applicable to any linear 
system, regardless of the degree of its characteristic polynomial, 
and is simple in application. It is worth while to note that the 
procedure presented here closely parallels the Nyquist procedure; 
however, whereas the Nyquist diagram presentation holds only 
for closed-loop systems, the present procedure was developed for 
open-loop arrangements. Finally, a word about construction of 
the frequency response curves; if the degree of the characterisic 
is an even number, the response curve will exhibit a cusp at the 
origin, while if the degree is an odd number, the curve will be 


smooth. 
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The Supersonic Thickness Drag of Rectangular 
Wings with Linear Spanwise Taper in 
Thickness 


K. Walker, Jr. 
Aerodynamicist, Douglas Aircraft Company, Inc., Santa Monica, 
alif 


March 30, 1953 
INTRODUCTION 


; SUPERSONIC THICKNESS DRAG for the rectangular wing 
with linearly varying thickness ratio (symmetrical about 
root chord) is presented in this note. The airfoil profiles con 
sidered have symmetrical modified double-wedge sections with 
constant thickness over a variable fraction of the chord length or 
biconvex sections. Recently, Henderson! made a similar calcu- 
lation for delta wings. 

Bonney? has calculated the drag coefficients for such rectangu- 
i.e., using the two-dimensional 
Strip 


lar wings using strip theory 
(infinite A.R.) drag coefficient for each wing section 
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theory is generally incorrect for those regions of the wing which 
are influenced by the root and tips; it is easily shown that strip 
theory is correct over the other regions of the wing. Strip 
theory results will be denoted by Cp. aaa 

In this note the exact corrections to strip theory are obtained 
using the well-known source-sink methods of linearized super 


sonic wing theory. 


(A) GENERAL MopIFIED DOUBLE WEDGE AIRFOIL 


The thickness at any point (y) along the mid-chord of a wing 
having a linear taper in thickness from root to tip may be ob 


2,6 ”] 


(1) 


tained from Fig. 1 as 


2, 0) te 


: oe 
¢(c/2, y) = ¢(c/2, b/2) + [(b/2) — y] 
“ . 5/2 


The slope of the surface of the generalized wedge may be calcu 


lated from Fig. 2, resulting in 


O2(x, y) ¢(c/2, y) ¢(c/2, b/2) 
Ox a! ye 


(b/2) — y ((c/2, 0) ¢(c/2, b/2) 
(2) 
b/2 ve 


Proper attention must be given to the sign of Eq. (2), since it is 
positive over the front wedge portion and negative over the rear 
portion. Care must also be exercised in the use of Eq. (2) in the 
region of the root Mach cone, since the thickness taper slope 
changes sign at the root. 

If the wing is extremely thin, it is permissible to evaluate the 
velocity potential in the plane of the wing ( = 0) for the pur 
pose of obtaining the pressure on the wing surface. This ve 


locity potential is 


(é, 7) 
: dit dn 


Oz 
uf f dé 
¢(x, y) = (3) 
T V (x — &)? By n)? 


After substituting Eq. (2) into Eq. (3), the potential at each 
point on the surface of the wing may be evaluated by integrating 
the source distribution over the fore Mach cone of the point. 
The drag of the wing may be written as 
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--=[" «fe » ax dy (4) 
Ox 


The derivatives of the velocity potential and Eq. (2) are now sub- 
stituted into Eq. (4) and integrated over the proper regions of 
the wing. Summing the integrated results gives the drag of a 
rectangular wing with a linear taper in thickness ratio and a gen- 
eralized wedge cross section. The equation for the drag coeffi- 
cient is 
7(3 — 4y) (1 - A) 

(A.R.pR)? 

471 — WF =) 

5 
m(A.R.p)? 


2 6, 
Cp = [iatataye? 
7 B 


where 6; = root thickness ratio, A.R.e = reduced aspect ratio 


(BA.R.), B = VM? — 1, and \ = ratio of tip to root thickness 
(6:/5;). The restrictions on Eq. (5) areO0 < y < '/s,2 < A.R.p, 
0 <A,andé,< 1. 


An equivalent thickness ratio may be defined by 
= 6,1 + A + A*)/3B = (5? + 5rd, + 5:?)/3 (6) 


where 6, and 6 are the root and tip thickness ratios. 
Eq. (6) to Eq. (5) results in 


2 8? 
*? 


Applying 


Cp = 


— 4y)(1 — A) 


| 4 ¥(3 _ 
(1 +A + A2)(A.R.R)? 


127(1 — y)%(1 — A)? 7) 
(1 +A + »*)e(A.R.p)? | * 


The drag coefficient obtained by strip theory is 


Co. a, = (2/7) (67/8) (8) 


Normalized drag coefficients will be found as the ratio of Eq. 
(7) to Eq. (8) in order to emphasize the errors involved in the use 


of strip theory. Two interesting cases are the double wedge 


(y = '/s) 
Cp 1—A 
—— 1+ e—») - 
Cre ar. 2(A.R.x)*(1 +A + A?) 
3(1 — A)? 
‘ —~ (9) 
27(A.R.2)9(1 + A + A?) 
and the modified double wedge (y = '/;) 
Ce ile 5(1 — X) 16(1 — 2)? 
Cc an WA.R.p)(L+A+A2)  9x(A.R.p)(1 +A +2?) 
(10) 


Eq. (9) has been plotted in Fig. 4 as the normalized drag co- 
efficient versus the spanwise thickness taper ratio for various re- 
duced aspect ratios. 


(B) Brconvex AIRFOIL 
The thickness at any point (y) along the mid-chord line may 
The slope on the surface of the bi- 
convex airfoil may be determined from Fig. 3 as 


be expressed by Eq. (1). 


Ox(x, y) _ ¢(c/2, v) I2 


Ox c/2 
2) (6/2) — y |] g(c/s =i 
. | b/2 l cl 


Following the same procedure 


— (4x/c)] = 


(11) 


as outlined for the generalized 
wedge profile, the drag of the rectangular wing with a linear taper 
in thickness ratio from root to tip and biconvex cross section can 
be written 


: 16 5,7] 1 (1 — ) 21 — rv)? ] 
Cp = - (1 +X 4A?) + 
, aE 5(A.R.p)? ~ 3(AR2) | 
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Restrictions implied on Eq. (12) are that 6, <<< 1,0 <A, and2 < 


A.R. pr. 
Using the definition for the equivalent thickness ratio defined 
in Eq. (6) and the drag coefficient obtained by strip theory 
Crow ax. = (16/3) (5/8) 
the normalized drag coefficient may be found as 
Cp 3(1 — A) 
=1+4+- ; = 
5(A.R. pz)? (1 + A + A?) 


2(1 — A)? 
m(A.R.p)(1 +A +? 
(13 


Coe A.R 


(C) ConcLupING REMARKS 


It is interesting to note from Fig. 4 that strip theory is adequate 
for engineering purposes provided the reduced aspect ratio is 
greater than 4. 


the pressure changes within the root and tip Mach cones cause 


For reduced aspect ratios of 2, on the other hand, 


the drag to be as much as 7 per cent higher than the strip-theory 
value. 
tained for all the profiles studied. 


Almost identical percentage drag corrections were ob- 
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The pressure changes within the tip Mach cone have the effect 
of reducing the drag below the strip-theory value, while the pres- 
sure changes within the root Mach cone have the opposite effect. 
The net result is to increase the drag above the strip-theory value, 
and the difference becomes larger with decrease in reduced aspect 
ratio. It should be understood that wings with inverse taper in 
thickness ratio (A > 1) have been excluded from the present study. 
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A Table of the Complete Cicala Function 


Yudell L. Luke and Dolores Ufford 
Midwest Research Institute, Kansas City, Mo. 
April 6, 1953 


_ AUTHORS HAVE SHOWN that the definite integrai 


F(x) = f ; e“(x +t —- Vx? + t*) dt/xt (1) 


can be completely expressed in terms of tabulated functions.! 
The integral was first introduced by Cicala? and has been called 
the complete Cicala Function. A short table was given by Cicala. 
Further entries were given by Jones.’ These entries are only to 
three decimals,and for the most part are recorded in Reissner‘ 
and Reissner and Stevens.. Aside from some unit rounding 
errors in all these tables, the entry for x = 0.05 in the latter two 
references is incorrect, but this does not affect the aerodynamic 
results. The integral appears in numerous studies dealing with 
the evaluation of spanwise loadings for oscillating airfoils.© 7 In 
view of the importance of this integral, recording of an extensive 
table permitting ready interpolation would be useful. 

We first bring together some relations describing the functions. 
Let the subscripts R and J denote the real and imaginary parts of 
Eq. (1). Employing the usual notation of Bessel functions (see 
reference 1), 

G(x) = Fr(x) — {1 — y — log 2x} = 


rf x 
- VJ, [Io(u) — Lo(u)| du + Li(x) — (x) ¢ (2) 


x 
F(x) = —(#/2) — (1/x) + Kil(x) + f Ko(u) du (8) 


where y is Euler’s constant. Power series expansions are readily 


obtained. Let 


A,~! = 22"+1 ym! (nm + 1)'(2n + 1) ! (4) 
B,-! = 2"+2(m + 1) (Qn + 3) {1 [nm +(3/2)]}?\ 
} 3 
C, = = + Dn 
2(n + 1) (2n + 1) - 
n (5) 
D, = > 1/m; Dy = 0 
m=1 
G(x) = (7/2) (= Anx**t! — , Byx?" *) (6) 
n=0 n=0 


x 


Fi(x) = (log 2 — y — log x) > 


n=0 


Anx™tl + 


x 


SY AnCax™*! — (4/2) (7) 


n=0 
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The function may also be defined by a differential equation. 


F(x) = if | T(u) — 1} (du/u?) (8) 
x 
T’(u) — T’(u)/u — T(u) = 1 u| 
j (9) 
7T(0)=1, T(0) = -i f 


x3F’"’ + 3x2F"” — x3F’ = 1 — ix (10) 


TABLE OF THE COMPLETE CICALA FUNCTION 


x G(x) Fr(x) — F(x) 
0 0 oo 1.5707 96 
0.01 0.0078 37 4.3426 45 1.5396 91 
0.02 0.0156 42 3.6573 02 1.5155 16 
0.03 0.0234 13 2596 O8 1.4939 56 
0.04 0.0311 51 9796 64 1.4740 94 
0.05 0.0388 7642 27 1.4554 938 
0.06 0.0465 31 5895 79 1.4378 96 
0.07 0.0541 72 .44380 69 1.4211 33 
l 
l 
1 
I 
] 
l 


ou 
“J 
tom wwe 


] 


0.08 0.0617 .3171 48 .4050 85 


o 2 
bo 
bo 


0.09 0.0693 59 2.2069 42 .38896 59 
0.10 0.0769 05 2.1091 28 .38747 87 
0.20 0.1506 66 1.4897 41 .2477 07 
0.30 0.2214 61 1.1550 71 .1461 74 


0.40 0.2894 53 0.9353 81 .0608 49 
0.50 0.3548 00 0.7775 84 0.9872 53 
0.60 0.4176 46 0.6581 09 0.9227 22 


.5644 40 0.8654 86 


bo 
a 


0.70 0.4781 


0.80 0.5363 75 0.4891 56 0.8142 80 
0.90 0.5925 02 0.4275 00 0.7681 56 
.00 0.6466 24 0.3762 61 0.7263 79 
.10 0.6988 46 0.3331 73 0.6883 63 
. 20 0.7492 65 0.2965 81 0.6536 30 
80 0.7979 74 0.2652 47 0.6217 86 
.40 0.8450 61 0.2382 26 0.5925 01 


.50 0.8906 05 0.2147 77 0.5654 93 
.9346 86 0.1943 20 0.5405 23 
70 0.9773 73 0.1763 82 0.5173 83 
.0187 35 0.1605 86 0.4958 93 
.0588 36 0.1466 19 0.4758 97 


ted) 


.0977 35 0.1342 25 0.4572 55 


ee ee ee) 


1 

1 

1 
2.10 1.1354 &9 0.1231 89 0.4898 45 
2.20 1.1721 52 0.1133 32 0.4235 59 
2.30 1.2077 74 0.1045 02 0.4083 01 
2.40 1.2424 02 0.0965 70 0.3939 84 
2.50 1.2760 80 0.0894 26 0.3805 32 
2.60 1.3088 53 0.0829 79 0.3678 74 
2.70 1.3407 57 0.0771 42 0.3559 48 
2.80 1.3718 31 0.0718 49 0.3446 98 
2.90 1.4021 14 0.0670 40 0.3340 72 
3.00 1.4316 34 0.0626 59 0.3240 25 
3.10 1.4604 27 0.0586 62 0.3145 14 
3.20 1.4885 22 0.0550 08 0.3055 00 
3.30 1.5159 47 0.0516 62 0.2969 49 
3.40 1.5427 29 0.0485 91 0.2888 29 
3.50 1.5688 95 0.0457 69 0.2811 09 
3.60 1.5944 68 0.0431 71 0.2737 64 
3.70 1.6194 72 0.0407 76 0.2667 69 
3.80 1.6489 30 0.0385 66 0.2601 O1 
3.90 1.6678 59 0.0365 20 0.2537 39 
4.00 1.6912 84 0.0346 27 0.2476 63 
4.20 1.7366 86 0.0312 39 0.2363 04 
4.40 1.7802 75 0.0283 08 0.2258 95 
4.60 1.8221 76 0.0257 57 0.2163 29 
4.80 1.8625 O8 0.0235 29 0.2075 13 
5.00 1.9013 76 0.0215 75 0.1993 64 
5.20 1.9388 70 0.0198 49 0.1918 14 
5.40 1.9750 83 0.0183 21 0.1848 02 
5.60 2.0100 90 0.0169 61 0.1782 73 
5.80 2.0439 68 0.0157 48 0.1721 81 
6.00 2.0767 82 0.0146 60 0.1664 85 
6.50 2.1545 52 0.0123 87 0.1537 47 
7.00 2.2268 80 0.0106 07 0.1428 03 
7.50 2.2944 54 0.0091 88 0.1333 03 
8.00 2.3578 43 0.0080 39 0.1249 84 
8.50 2.4175 23 0.0070 94 0.1176 38 
9.00 2.4738 95 0.0063 08 0.1111 06 
9.50 2.5273 02 0.0056 48 0.1052 60 
10.00 2.5780 32 0.0050 84 0.0999 98 
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1 — ixF(x), where 


An alternative definition is S(x) = 


S(x) 


= (inx 2) fo [H\(iu) — Y\(in) — (2/r)]| (du/u) 
x 


The asymptotic expansions take the form 


Fp(x) ~ [((1/2)x?] D5 x-™ [1-3-5 .. 


n=0 


.(2n — 1)]? X 


+ 1)/(n + 1) 


(2n 


F(x) ~ —(1/x) + (9x/2)'/? e-* [1 — (9/8)x 4 


(345/128)x? — |/x? 


The elements in the expansion of Eq. (13) satisfy the recurrence 


relation 


16(m + 2)dm42 + 6(2m + 3) (2m + 5)dm4i + 


(2m + 1) (2m + 3) (2m + 5)am = 0 (14) 


Tabulated here are values of G(x), Fr(x), and F/(x) for x = 
0(0.01)0.1(0.1)4.0(0.2)6.0(0.5)10.0 to six decimals. The 
tions were calculated using the tables cited in reference 1 where 
applicable. 
solution of pertinent differential equations. 
of a table of 


func- 


Otherwise, computation was by series or numerical 
Use was also made 


_ Px 

Ko(x) = f, Ko(u) du 
constructed by the authors.’ The errors in the G(x) and F/(x) 
tables do not exceed 2 and 0.6 units in the last decimal place, re 
spectively. 

Though G(x) is given for the entire range, it was particularly 
designed to facilitate interpolation of F(x) in the neighborhood 
When interpolating in the G(x) tables, virtually 
full accuracy is assured if five-point Lagrangians are employed. 
The same is true for Fr(x) and F(x) if x > 1.0. 
the latter in the vicinity of the origin is facilitated using the table 


of the origin. 
Interpolation of 


of Ko(x) mentioned above. In this range, the series expansions 
Four terms of the asymptotic expan- 


This 


are rapidly convergent. 
sion for F(x) gives at least five decimals accuracy if x > 7. 
is also true for F(x) if Eq. (13) is employed as written. 
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SCIENCES 


The tables of F(x) both overlap and supplement the values 
given by Dingle and Kuessner,’? who tabulated S(x) to five deg. 
mals. Comparable values agree. These authors give ay 
asymptotic expansion equivalent to Eqs. (12) and (13), but for 
Eq. (13) only the term 1/x appears. 

In some recent studies,® 7 need has arisen for values of the 
incomplete Cicala Function, and some values of this expression 
have been recorded. It can be shown that tabulation of the 


incomplete Cicala Function can be reduced to tabulation of the 


y , 
f, (e“ — 1) Vu? 
0 


Tables of the latter'® are available to seven decimals for 6 = 
0(0.2)1.8, y = 0(0.05)2.0(0.1)2.5; b = 2.0(0.2)4.0, y = 0(0.05) xX 
0.5(0.1)2.0. 


functions 


+ b? (du/u) 
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